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Preface 



THESE LECTURE NOTES are the outcome of a six week course with 
the Bangalore applied mathematics group of the Tata Institute. It is a 
pleasure to recall the enthusiasm and energy of the participants. The 
lectures were on introduction to certain aspects of gas dynamics con- 
centrating on some of the currently most important nonlinear problems, 
important not only from the engineering or computational point of view 
but also because they offer great mathematical challenges. The notes, I 
hope, touch on both these aspects. 

I am indebted to Professor K.G. Ramanathan for inviting me in the 
first place and for making my visit so enjoyable and stimulating. To 
RS. Datti, who cooperated in writing these notes and worked and re- 
worked them, go my very deep thanks. I would also like to thank RR 
Gopalakrishnan who reproduced them so well. 

Cathleen S. Morawetz 

New York 
December 1980 
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Chapter 1 



The Traffic problem and a 
first order nonlinear 
equation 

1.1 Introduction 

Many physical laws occur as conservation laws. The most general such 1 
law in its differentiated form is given by 

Y, + (F(Y)) X = 

where Y - (Y\, . . . , Y n ) is a vector valued function of x, t with x e R" 
and F(Y) is a matrix valued function of Y. The term conservation comes 
from the fact that if F(Y) — > as \x\ — > oo then J F|<ix| is constant for 
all time, that is, thes integrals are conserved. 

In these notes we shall consider conservation laws in only two in- 
dependent variables, either time and one space variable or two space 
variables. 

The simplest case to treat is that of the single conservation equation 
with the variables t, x. Then Y, F(Y) are scalar. We replace Yby N and 
consider 

N t + (F(N)) X - 0, -oo<x<™, t>0 (1.1) 
1 



2 
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with N(x, t) e R. 

We shall now formulate the traffic problem, first proposed, by 
Lighthill and Whitham [26]]. Let N{x,t) denote the density, the num- 
ber of vehicles passing through the position x at time t on a highway. 
Let u(x, t) be the average (local) velocity of the vehicles. Then in any 
section the conservation equation states that the total number of 

cars is preserved, or 



2 Assuming the quantities N, u to be smooth, we obtain in the limit t\ — > 
t 2 , that 



This is the integrated form of the conservation equation. As x\ — » x%, 
we obtain 



We can always put this equation in the form of (11.lt if we use u - U(N). 
This assumption seems to be reasonable since drivers are supposed to 
increase or decrease their speed as the density N decreases or increases 
respectively. The maximum value of u occurs when N = (the maxi- 
mum is, say, the maximum allowed speed). And when N is maximum 
u = 0. Hence the graph of U plotted against ,/V takes the form shown in 
figure 1.1. 






N, + (Nu) x = 0. 



(1.2) 
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Rewriting (11.2ft as 



N t + (F(N)) X = 0, 



where F(iV) - NU(N) is the flux of cars, we see that F - when N - 
and when jV is maximum (in this case U(N) = 0). Hence the graph of 3 
the flux curve F(N) will look as in figure 1.2. It could have the shapes as 
shown in figure 1.2 (a) and 1.2 (b), but these lead to certain difficulties 
in the theory. 




Fig. 1.2. 



4 



1. The Traffic problem and a first order nonlinear equation 




Fig. 1.2 (a). Fig. 1.2 (b). 

Consider the simplest case in which the flow is steady, i.e., indepen- 
dent of time. Then from (1.3), we obtain 

F(N) - constant. 

In general, the line F(N) - constant, will cut the flux graph in two points 
as shown in figure (1.3). 




Fig. 1.3. 



4 If we require continuity in N, then we have to take either N - Ni or 

N = N2- If we allow jump discontinuities in N, then uniqueness in the 
solution will be lost as shown in figure (1.4). Then thick line can be a 
condidate for the solution but the dotted line could also be a candidate 
for the solution. 

So, for uniqueness, we need an additional condition which is called 
an Entropy Condition. The terminology will become clear when we 
study gas dynamics. Again this has to come from the physics of the 
problem at hand. 



1.1. Introduction 
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For the traffic problem at hand, we would like to add the entropy 
condition that infinite acceleration is impossible, i.e. 

dU , ON 

— < oo(or — > -oo). 

ox ox 

It turns out that there is then a unique solution to the initial value prob- 
lem. 



) l N 



No 



iVi 



Fig. 1.4. 



We also observe that for a given conservation law in differentiated 
form there are several equivalent conservation laws in differentiated 5 
form. For example, consider 

N t + (F(N)) X = 0. 

Let P(N) be any arbitrary integrable function of N. Put 

N 

Q{N) = J P{M)dM. 

N, 

Then 



Q: = Q'(N> -N, 
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N 

= ~(R(N)) X , where R(N) = J P(M) ■ F\M)dM. 

AT, 

Thus 

Qt + R x = 0. 

Hence to choose the correct entropy condition, and thus to get uni- 
queness, we should look at the integrated form of the conservation law. 

If we allow discontinuities in the solution, it is a so-called weak so- 
lution. We shall now give a precise definition of a weak solution for a 
general first order nonlinear equation and then return to the traffic prob- 
lem again. 



1.2 Weak Solutions 

Consider a first order nonlinear equation 

N t + (F(N)) X - 0, -oo<;t<oo, t>0. (1.4) 

Let Af be a classical solution of dl.41 ). By this we mean that N is a C 1 
function in x, t variables and satisfies (11.4ft identically. Let now x(x, t) 
be any C^-function (even C 2 is enough) which vanishes for large \x\ and 
t and t = 0. 



X(x, t){N t + (F(N)) x }dt dx = 0. 
Integrating by parts, we obtain 

(X,N+XxF(N))dtdx = (1.5) 



since the boundary terms vanish. Also, the entropy condition — > -oo 

ox 

becomes 

oo oo 

(f)(x, t) ■ —dt dx > -oo 
ox 



oo oo 
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for all cD > o, C°° function vanishing for large \x\ and t and t - 0. This 
again by integrating by parts, becomes 



OO OO 



$> X N dt dx < oo. (1.6) 



Motivated by this, we now define a weak solution. 

Definition . A locally integrable function N(x, t) is a weak solution of 
dl.4t if dl-5l > and ( 11.61 ) hold under the conditions stated there. 

Note that a classical solution (usually called strong solution) is nec- 
essarily a weak solution. Conversely, it can be seen that if ,/V is a weak 
solution and N is of class C 1 then N is necesaarily a strong solution. 



1.3 Initial value problem 

We now consider the problem of solving a general first order nonlinear 
equation 

N, + (F(N)) X = 0, -oo < x < oo, t > (1.7) 
with initial data 7 

N(x, 0) = cD(x), -oo < x < oo. (1.8) 

If we put G(N) = F'(N), < fT77l can be written as 

N, + G(N)N X - 0. (1.9) 

This we solve by the method of characteristics. If we define a curve C 
in x, t plane by dx/dt = G(N), we find that on C, (11.91 reduces to 

dt 

or N = constant on dx/dt = G(N). Since this means G(N) is also a 
constant we see that the characteristics are straight lines with slopes 
(G(AO)" 1 . If x(0) - £ we find that 

N(x, t) = <&(£) on dx/dt = G(N) 
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with x(0) = 

Hence we have a solution in implicit form given by 

x = £ + fG(<D(£)) 
N(x,t) = ®(0 

If we can determine £ = £(x, t) from the first equation, then we know N 
at (x, t) uniquely. But, however, this is not always the case; and trouble 
occurs if the characteristics intersect, as we shall see now. If £ = R, 
t; - L are two points on the x-axis such that 

G(<t(L)) > G(<D(/?)) > 




Fig. 1.5. 

(See figure (1.5)) then the characteristic through R intersects the 
characteristic through L at P. Thus the value of £ corresponding to the 
point P is not unique. 

However, if the characteristics fan out, then clearly there will be a 
unique characteristic through every point and the solution will be deter- 
mined uniquely. If 

G(0(L)) < G(0(/?)) (1.11) 

for L < R then the characteristics in the x, t plane have decreasing slopes 
and they never intersect. In such a case, we obtain a unique continuous 
solution. 

Since the solution is given implicitly by d 1 - 1 Oft another way of look- 
ing at the breakdown, when dl.llt does not hold is to try to solve the 



(1.10) 



1.3. Initial value problem 
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first equation in dl.lOt using implicit function theorem, the functional 
relation 

/(*,*,£ = 
can be represented in a single- valued way by 

if and only if / f 56 0. Put 

f(x,t,Q = tGm&)-{x-&. 



Then 



Hence, if 



dG dO , 

ff-t — + 1. 

H dN <% 



dG dO 
dN ' ~dj 

is always positive for all positive t we have fg > 0. On the other hand, if 
this expression changes sign, there is a finite time, called breaking time, 
at which fg = o. We note from the relations 

AJ _ <&' ■ G(N) M G(N) 



1 + tO'G' 1 + tO'G' 

that the derivatives in ,/V will blow up at the time of breaking. We also 
note that the breakdown in the solution can happen even if the initial data 
are very smooth. Suppose the flux curve is convex, so that G' = F" < 0. 
Now, if the initial data O is very smooth and tends to zero as \x\ — > 00, 
then will change its sign and there is always a time at which the 
solution will become singular. 

This is true in higher dimensions also. John, F. [19] has investigated 
the question of how late with 3-space variables and a nonlinear equation, 
a breakdown can occur. Note in the above that as |0'| — > the breaking 
time t — > 00. Klainerman [21] has shown for the nonlinear analogue of 
the wave equation. 

d 2 N 

ij 1 7 
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that if the number of space variables is greater than or equal to 6, then 
compact initial data can propagate smoothing for all time. This leaves 
open the anologous important question in 2 and 3 dimensions. 

Exercise 1.1. Describe the solution of 

N t + (F(N)) X = 

with 

F(N) = 16 -N if N<4 
= if N>4 

for the initial data 

a) N(x, 0) = 2 - tanh x 

b) N(x, 0) = 2 - tanh x if x < 

= 2 tanh x if x > 0. 

2 

Where does the discontinuity in N x move ? 

Exercise 1.2. Work out the corresponding theory for 

U t + a(U)U x + b(U)U y = o 

and determine conditions on a and b which lead to discontinuities for all 
compact initial data. 

We now return to our traffic problem and ask how the condition 
dl.lU should be interpreted as a condition on the density at some given 
time. 

We refer to the flux curve (see figure (1.2)). Suppose F attains a 
maximum value at density Af*. Then F is increasing in [0, ./V*] and de- 
creasing in [N*, AWx]- A^Max is the maximum density at which F = 0. 
Hence G(N) = dF/dN is decreasing in [0, N t ], equal to sero at and 
again decreasing (negatively) beyond A^* . Hence the graph of G(N) will 
look as shown in figure (1.6). 



1.3. Initial value problem 



11 




If the traffic problem is accelerating at the given time, i.e. cars to the 
right going faster than those to the left then N is a decreasing function 
of x. Thus 

N(L) > N(R) 

for all pairs L, R such that L < R and thus 

G(N(L)) < G(N(R)) 

(unless the density has increased beyond N t ), so that fll.lU holds and 
we have a unique continuous solution. However, the tail of the G(N) 
curve for N > N* was fixed up artificially to give us a smooth curve and 
we should properly speaking ignore this region. 

If the traffic is decelerating then N(L) < N(R) for all pairs L, R with 
L < R and we will not be able to obtain a continuous solution. 

We are then forced back to re-examine our model. The conservation 12 
of the number of vehicles still holds, but we can allow discontinuities in 
density. 
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We have the conservation law 



(Ndx - Fdt) = 



e 



for every closed curve e in the x, t plane or equivalently 




N(x, t)dx + [F(N)]* 2 = 



(1.12) 



for every segment [xi,X2] if N has jump discontinuities. We ask for 
piecewise smooth solutions N and investigate what happens across a 
discontinuity in ,/V as in the steady case. Such a discontinuity is called a 
shock and the curve on which it lies is the shock curve. 

A discontinuity in N implies, of course, a discontinuity in the speed 
U. These discontinuities represent instant or rather infinite decelera- 
tions. In many situations, the presence of a discontinuity represents a 
pile-up and to say the least, the model breaks down. However, densities 
have been defined in an average sense and if the traffic is this enough 
a sharp deceleration can occur and we can study what happens to it. 
Ofcourse, the relation of speed to density inside this narrow region rep- 
resented by the discontinuity cannot be the old one. 

On the other hand, one might also ask, why not allow discontinuities 
in an accelerating situation ? It is not quite clear that in a traffic situation 
13 one should not but again it involves violating the speed density relation 
in some narrow region and this time we would lose uniqueness. 

Let x = s(t) be a C 1 -function representing the shocks curve. We 
work with the conservation law given by dl.121 ). Let x\ < s(t) < X2 at 
some time. Then from dl.l2t we obtain 



s(t)- 





Ndx + 



s(t)+ 



Differentiating under the integral sign and taking the limits 



xi — > s(t)— 
x 2 -* s(t)+ 



1.3. Initial value problem 
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we obtain 

S(N r -Nl) + (F r -F l ) = 0, 

where we denote ds/dt by S and Nr, Nl and Fr, Fl are the correspond- 
ing limiting values from right and left respectively. We then have 

_ F R -F L 
Nr-N l 

S is called shock speed and (Nr - N£) shock strength. The equation 
d 1 - 1 3I > is called a jump condition. 

Exercise 1.3. If ./V is a piecewise smooth weak solution of N t +(F(N)) X = 
0, show that the jump condition holds across the line of discontinuity. 
Further, show that as the shock strength goes to zero, the shock speed 
becomes the characteristic speed. 

We want to show that by allowing shocks we can solve any initial 14 
value problem uniquely and that we can study in particular where the 
shock travels and how strong it gets. We first look at constant states. 

The simplest problem to solve is the transition from one constant 
density, say, Nr to another constant density Nl- The deceleration re- 
quirement is that Nr > Nl- The shock speed, by (11.131 . is the slope 
of the line segment connecting two points on the flux curve; see figure 
(1.7). 




*~N 



Fig. 1.7. 



It is positive or negative depending on the relationship of Nl and Nr 
to N*. A more important question is whether it is more or less than the 
speed of the traffic. 
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Generally speaking, relative to the traffic ahead at time t, the shock 
is always retreating, i.e., S < Ur. For, since F - NU. 

S ~ U R = 7F " T } ~ u * 
N r -N l 

N l (Ur ~ UP _ . M M 

= ~^Tr ^7 < °> smce N * > N l- 

Nr-N l 

15 Similarly, 5 < Ul- Thus, all traffic ahead of the shock remains ahead but 
all traffic behind it eventually hits it and decelerates. The path history of 
cars is illustrated in figure (1.8). 




1.4 Initial value problem with shock 

We already know from the constant configuration (and it can be proved 
generally) that the vehicles after crossing the shock move off at a slower 
speed. However, the behaviour of characteristics is quite different; they 
hit the shock from both sides. In the x, t plane the slope of the charac- 
teristic satisfies dx/dt = G(N) = dF/dN. From the G - N curve (figure 
(1.7)), we see that not only G(N R ) < 0, but G(N R ) < S , for N > N* since 
Nr > Nl and F is convex. Hence the characteristics ahead of a shock 



1.4. Initial value problem with shock 
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run into it. Behind the shock the reverse is true, G{Nj) > S and again 
the characteristics hit the shock (figure (1.9)). 



t 



Shock 




Fig. 1.9. 



If we knew how to lay down the shock across overlapping character- 
istics, we could solve the initial value problem. The shock starts when 
the continuous method breaks down, i.e., when dN/dx becomes infinite 
at some time and plane. Let this, for the sake of argument, be (0,0). 
Then initial slope of the shock is the characteristic slope and we inte- 
grate from this point (0,0). The shock velocity 

dx _ _ F R - F L 
dt~ " N R -N L 

where CD is the given initial data and %R(x,i),{;L,(x,t) are the abscissae 
of the appropriate right and left characteristics. To find this curve, a 
useful approximation is often made. We use the subscripts R and L in an 
obvious manner and note two possibilities of getting Taylor expansions 
for the shock speed in terms of the shock strength Nr - Nl- 

F — F dF 1 d^F 
S = A^vf = + 2 { df )R{NR ~ NL) + ° ({Nr ~ NL)2) 

= ( J^k + \(^)l(N l ~ N R ) + 0n((N L - Nr) 2 ) 



16 
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as Nr — > Nl- Adding the two expressions, 



dNj R \dNjJ 2\\dN 2 ) R \dN 2 
(N R - N L ) + 0((N R - Nl) 2 ). 



d2p \ 



or, expanding ^^jj an d yjpj2 j * n terms °f (Nr - Nl), we see that 



>L 

This means the shock speed is the average of the slopes of the two char- 
acteristics with a second order error. Note, if F is quadratic the formula 
is exact. 

To find the whole flow in the (x, 0-plane, we now have a relatively 
simple reciple: 

1. Carry Af along the characteristics coming from the initial line, i.e., 
lines of slope (dF/dN)~ l - G~ l . 

2. Note the region where they are crossed. 

3. From the points on each line t = constant where such a region 
starts, start a shock. 

4. Integrate the differential equation for the shocks, 

dx 

-dt =8{XjX 

where g(x, t) = -(Gr + Gl) is a good approximation. 

This procedure is quite straight forward until shocks intersect. But 
there the initial value problem can be restarted provided we can solve the 
local problem with constant states on each side. Note if initial data have 
a discontinuity then the two sides can be connected either by a shock 
or by a degenerate solution of (11.101) with x - tG(N) called a centered 
wave. However, interaction shows that the number of shocks decreases. 



1.4. Initial value problem with shock 
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The real numerical difficulty of the foregoing procedure lies in in- 
verting the variable and finding £ as a function of x, t. For this reason, 
it is often preferable to use a difference scheme. Before going into this, 
we now give some existence and uniqueness theorems. 

Remarks . We note the difference in shock speed if we use different 
conservation laws that are equivalent in the smooth case. From 

N t + F x = 0, 

we have seen that, Q satisfies 

Qt + T x = 

where 

N 

Q = J P(M)dM 

and 

N 

T = J P(M)F'(M)dM. 
Considering shock speeds, we see that 

s = Tr- T l 
Qr-Ql 

Nr 

J P(N)F'(N)dN 

_ Nl 

J P(N)dN 

N L 

In the limits as Nl — > Nr we see that 5 tends to the characteristic speed, 19 
the same for both the forms. But for Nr + Nl we can obtain a wide 
range of shock speeds by the choice of P(N). 



To get uniqueness, as we shall see, we introduce the following en- 
tropy condition to go with a specific weak conservation law: 
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The characteristics starting on either side of the discontinuity curve 
(shock curve) when continued in the direction of increasing time t inter- 
sect the line of discontinuity. 

This means the characteristics issuing from a shock go backwards 
in time (see figure (1.9)). This will be the case if 

F\N L ) > S > F'(N R ) 

where Nl, Nr are the values of N on left and right of the shock curve 
respectively and S is the shock speed. 

Hereafter, in this section, by a shock, we mean a discontinuity satis- 
fying the jump condition and this entropy condition. For the following 
results, which we are going to derive, we refer to P.D. Lax E31 . The 
most important feature there is the role of the L 1 norm. A piecewise 
continuous function such as a solution with shocks will have its first 
derivatives in L 1 but not in L 2 . Consider 

N t + F x = 0, -oo<x<oo, t>0. (1.14) 

We assume F is convejfl and a C 2 -function. 

20 Exercise 1.4. Let / : (a, b) — > R be a C 2 -function. Then / is convex iff 
f"(x) > for all x e (a, b). Further show that / satisfies the inequality 

f(x)>f(y) + (x-y)f'(y) (1.15) 

for all x,y € (a, b). 

Theorem . Let N, M be two piecewise continuous solutions of M.14\ 
whose discontinuities are only shocks. Then \\N — M\\(t) is a decreasing 
function oft where 

CO 

\\N-M\\(t) = J \N(x,t) - M(x,t)\dx. 

— oo 

(We assume the integral exists). 



'In the traffic problem F is concave but the problem can be transformed to this case. 



1.4. Initial value problem with shock 
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Corollary (Uniqueness Theorem) IfN = M at time t - 0, then N = M. 

Proof. Let x - y n (t) be the points such that (N - M) has the sign of 
(-1)" my n {t) <x< y n+ i(t). Then 



y n +\ 

l.x. 

yn 

There are two cases: 



\\N - M\\{t) = f(N- M)di 



(i) Suppose N = M on a curve which is not a shock curve of either 
solution. Then G(A0 = G(M). Thus the two sets of characteristics 
have the same slopes and hence coincide. Hence there is a seg- 
ment or a point on every line t = constant where N = M. If it is a 
point, the curve on which N = M is a characteristic; if it is a line 
segment then the whole region swept out by the characteristics 
from the segment satisfies N - M. 



(ii) The curve where N = M is a shock. Consider 21 

t )dx+ 



j t \\N-M\\(f) = Y J (-lf 



yn 



Consider the term in the bracket; from equation d!.14t it is 

i 

(F(M) x -F(N) x )dx+(N-M)^-%: +l - [ F (m)-F(N)+(N-M)^-] 



y n +i 

v „ , x , v , , x ,^ , v , ^ , ^ lyn - l x v»v - v w v - dt iy n 

y n 

(1.16) 

where y n are now points of discontinuity of Af or M. The contribution 
from other points y n vanishes because N = M and F(M) = F(N). Next 
we calculate the contribution at the upper point y„+i. Suppose has 
a discontinuity at y n+ \ and M does not; the other cases are analogous. 
Now 

dy = F(N L ) - F(N R ) 
dt N L -N R 
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Since (N - M) changes sign and M does not have a discontinuity at y n+ i, 
we must have 

N R < M < N L . 

The contribution of the term in the bracket in d!.16t at y n+ \ is thus, with 
N = Nl, given by 

F(N L ) - F(N R ) 



F(M) - F(N L ) + (N L - M) ■ 



N L -N R 



{ M-N R n l -m 

< 0, (by the convexity of F since N R < M < N L ). Since M < N L , 
(Nl - M) is positive and therefore n is even and hence the contribution is 
negative. Arguing on the same lines, we find similarly that the contribu- 
tion from the lower point y n is also negative as well as the contribution 
when both Af and M have shocks. This completes the proof. □ 

Remark. A similar estimate which yields uniqueness can be made un- 
der alternative conditions on F, other than convexity. 

Theorem (A minimum principle) Consider the initial value problem 

N t + F x = 0, -oo < x < oo, t > 0, 
iV(jt,0) - <D(x). 

Let N(x,t) be a continuous and differentiable solution. Let <D € L' (It 
suffices to assume that O vanishes for large negative x). Put 



I(x,y, t) = I ®(s)ds + tH(^—^), 



where 

H(L) - MG(M) - F(M), M - G~ l (L), G = dF/dM. (1.17) 



, x — y 

Then N(x,t) = G ( — - — ) where y minimizes I(x,y,t). (Note G(N) - 

F'(N) is an increasing function of N by our assumption on F; so that 
the inverse G~ l exists). 
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X 

U(x,i) = J N(y,t)dy, (1.18) 

— oo 

then 

U X = N. (1.19) 
Since N satisfies the differential equation, we obtain 

U, + F(U X ) = 0; (1.20) 

here we have adjusted the integration constant by putting 

F(0) - 0. (1.21) 



If (I1.15t is applied with U x and any number M, we obtain 
F(U X ) > F(M) + (U x - M)G(M). 

Or using dl.20t . 

U, + G(M)U X < MG(M) - F(M). (1.22) 

Let y denote the intercept on the x-axis of the line given by dx/dt = 
G(M) or 

{x -y)/t = G(M), t>0. (1.23) 
Integrating ( 11.221 1 along this line w.r.t. t from t - 0, we obtain 

U(x, t) - U(y, 0) < t[MG(M) - F(M)] . ( 1 .24) 

From (11.231 . we have 

G _1 (^p) - M. (1.25) 



If H is defined by d 1 - 171) . we see that 

U(x,t) < U(y,0) + tH{^—^) = I(x,y,t). (1.26) 
We also note that 24 
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dH/dL = G~\L). 

Let G(0) = c; then G~\c) = 0. Since F(0) = 0, we have H{c) = and 
this is its minimum value. The inequality in (ll.26t holds for all choices 
of y. In particular for the value of y for which M, given by (ll.25t . equals 
N(x, t), the sign of equality holds in (11.241) along the whole characteristic 
dx/dt = G(N) issuing from (x, i). Therefore, the sign of equality holds 
in (U.26t . This completes the proof. □ 

Remark 1. The above theorem holds also for generalised (weak) solu- 
tions whose discontinuities are shocks. For, relation ( 11.201 ) is the integral 
form of the conservation law and so relation (11.261) is also valid for gen- 
eralised solutions. Since all discontinuities are shocks every point (x, t) 
can be connected to a point y on the initial line by a backward char- 
acteristic (entropy condition). For this choice of M equality holds in 



Remark 2. The converse of the result given in remarkfTJis also true. 



An estimate for large t. In the first theorem, we have found an explicit 
formula for the solution of an initial value problem in terms of its initial 
value. Recall _ 

N(x,t) = G~ 1 (^) (1.27) 

25 where y minimizes 

y 

l{x,y,t) = J ®(s)ds + tH(?-y-). (1.28) 

— oo 

We have also seen that H takes its minimum values at F(0) = c amd 
H(c) - 0. Let 

k - \{G-')\c) = \h"{c). (1.29) 

Assuming F is strictly convex, we have k > 0. Suppose there are two 
positive constants k\, &2 such that 



2*i < H" < 2k 2 . 



(1.30) 
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It follows then from H(c) = 0, H'(0) = that 

H(L) > k x {L-cf. 



Thus 



tH( -) > — (x-y-ctf. 

t t 



(1.31) 



Now we assume $ e L 1 and let ||<1)|| denotes its L'-norm. Then using 
( 11.28b and dl.311 . we conclude 



-\m\ + j(x-y-ct) 2 <I(x,y,t). 



(1.32) 



But 



x—ct 

I(x, x-ct,t)= J $>(s)ds < ||0|| 

— oo 

and hence by the minimum principle 

/(*,?,*) <||<D|| 

at the minimizing function y. Combined with dl.311 . we obtain 

1/2 



x-y /2||0||V /Z r 

c\ < < > -m/yt, say. (1.33) 

t { k\t ) 



From (G -1 )' - H" < 2k 2 and G~\c) = 0, we have 
\G~\L)\<2k 2 \L-c\, 



or 



\G-\^)\< 2k 2 \^ -c\ 



Thus 



< 2k 2 ■ ml Vf = mj Vf, say. 



\N(x, t)\ < ml V?, from Q7b 
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Now suppose that O vanishes outside (-A, A); then 



; 



<£>(s)ds - 0, if x, -A 



= constant, if x > A. 
According to (11.33b . the minimum value of y lies in the interval 
x-ct-m^Jt<y<x-ct + m Vf- 

If x < ct - m y/i - A, then y < -A, where the value of 

y 

J ®(s)ds 

— oo 

is independent of y and therefore the minimum of / is attained at the 

x — y 

point which minimizes tH{ — - — ); this value is y - x - ct. Similarly, if 

x > ct + m yjt + A, the minimizing value is y = x - ct. Since G _1 (c) - 0, 
we conclude from N(x, t) = G _1 (^r") that N(x, t) = for x outside the 
27 interval (ct — msfi — A, ct + fhyft + A). This can be restated as: Every 
solution ,/V whose initial value vanishes outside a finite interval, at time t 
vanishes outside on terval whose length is 0( V0 and inside this interval 
it is 0(1/ V?). 

In fact, the result can be improved. We state the theorem without 
proof. 

Theorem. Define the 2-parcimeter family of functions M(p,q), p,q > 
as 



M(x,t; p,q) 



\{x - ct)/tF"(0), for ct - ^Jpi < x < ct + ^fqt, 
0, otherwise 



Let N{x, t) be any solution with shocks of 

N t + (F(N))x - 
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where F is convex, F'(0) - c. Then 

\\N-M(p,q)\\(t) 

tends to zero as t — > oo where \\ ■ \\ is the L l -norm and 

y 

p = -IF" '(O)min J ®(s)ds 

— CO 
CO 

4 = -2F"(0)max J ®(s)ds. 

y 

1.5 Modifications by diffustion and dissipation 

Consider the first order nonlinear equation 

N, + NN X - 0, -oo < x < oo, t > 0, (1.34) 
JV(;c,0) = <D(jc), 

whose solution is implicitly given by 

* = f + rt>(^ 
N(x,t) = <s*g). 

We have seen that however smooth CD may be the solution becomes 28 
discontinuous after a finite time, if O has compact support. We want to 
consider a different model in which we add an extra dissipative term: 

N t + NN x = aN xx , (1.35) 

where a is a small positive number. This is the well-known Burgers' 
equation that has been utilized by Burgers as a mathematical model of 
turbulence. By a nonlinear transformation (Cole-Hopf transformation), 
it is known that ( 11.351 ) can be reduced to the linear heat equation 

% - (1-36) 
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where 

vp 

N = -2a—. 

Solving (11.361) a complete solution of dl.35t with N(x, 0) = O(x) is given 
by 

CO 

/ ^e-x /2a dn 
J e-x' 2a dn 

— oo 

where 

n 

Xim x,t) = J <btf)dif + (x - nf/lt 

o 

Hence in the case of Burgers' equation a solution exists for all time ever 
29 for discontinuous (bounded, measurable) functions O. It can be shown, 
by the method of steepest descent, that as a — > the solution of the 
Burgers' equation becomes, asymptotically, a solution of (II .34ft . 

Now instead of putting an extra dissipative term, we put a dispersive 
term and consider 

N t + N xx =pN xxx . (1.37) 

This is the Kortaweg-deVries equation, originally formulated in the 
study of shallow water theory. Recently, P.D. Lax and D. Levermore 
have developed the theory of the connection between this equation and 
shock theory in the limit f5 — > 0. Many other mathematical properties of 
this equation have been studied in recent years. 

1.6 Propagation of singularities in derivatives 

Let Q. = Q.i be an open set in (x,t) plane (see figure 1.10). Let 
f(x,t) be a continuous function in Q.. Suppose f x , f t are continuous in 
Oi and Q2 and have finite limits as they approach T from either side. 
Let T be described by a smooth curve 
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1 i 




Fig. 1.10. 



x - sit). Let [•] denote the jump across T. Then a calculation shows 30 

4 LA = [/^' + L/il- (1-38) 

at 

Consider now the equation 

N t + G(N)N x = (1.39) 

in Q. Let N be a continuous function in O, having discontinuities in N x , 
N t across T, and satisfying dl.39t in Q.\, Q.2- Applying d 1 .38l > to N, we 
abtain 

[N x ]s' + [N t ] = 0, (1.40) 

since [N] = 0. 

Since N satisfies dl.39l > in Q.\ and Q.2, we obtain as we approach T 

that 

[N t ] + G(N)[N X ] = 0. (1.41) 
Let A = [N x ]; then from (fL4~0l [JV t ] - -As*. Thus from OB, we obtain 

-As' + G(A^)/l - 0. 

Assuming A + 0, we obtain 5' = G(N). But then this curve is precisely 
a characteristic. Hence the discontinuities in N x , N t propagate along 
characteristics. 

Note. This is true in higher dimensions also. 
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Let us look more closely at the singularities in N x , N t . To do this, 
we put p = N x , q = N t . Then in the regions where p, q are continuous 
they satisfy the equations 

Pt + G{N) Px + G'{N)p 2 = 0, (1.42) 
q t + G{N)q x + G'(N)pq = (1.43) 

31 respectively. Solving (11.421) by the method of characteristics, we find 

p = on dx/dt - G(N), 

J G'(N)dt - c 

where c is a constant. Note that the characteristics for these equations 
and the original equation are the same. Note that since at the time of 
breaking p becomes infinite, the constant c must become the integral in 
the denominator at that time. 



1.7 Computing methods 

Although we have obtained a method of solving a nonlinear equation, 
it may be difficult to obtain the solution explicitly using the method 
described. In practice, it is preferable to use a difference scheme. 

We first consider a linear equation 

U, + aU x = 

where 'a' is a constant, with initial condition U(x, 0) = $>(x). Let the 
domain be approximated by a rectangle. Divide this rectangle into small 
rectangles of length '/z' and width (see the figure 1.11). We want to 
find an approximate value of U on this mesh of points. 
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Fig. 1.11. 

Let U(x,t) - U(hi,kj) - Uj. The simplest way to approximate 
dll /dt, at (hi, k j), would be by 



U 



7+1 



U J 
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making an error of order k. We shall this is too inaccurate. If a similar 
approximation is made for dU/dx we obtain the following difference 
equation 



uj +1 - uj 



+ a 



U J ., - U{ 

i+\ h 



-0 



or 



U; 



i+i _ 



ak 



= u> - -(u J M - Uj) 



(1.44) 



Its simplicity lies in the fact that it requires the values only on the first 
row, which will be given by the initial data. Each evaluation, however, 
gives an error of order h 2 or k 2 - a 2 h 2 , say, and if one substitutes con- 
secutively, the value U. is obtained and involves using the approxima- 
tion j(j - l)/2 times, i.e., making an error of order j(j - l)h /2. Now, 
if j is a mesh corresponding to final time T, j = T/k = T/h and, hence, 
the error behaves like T 2 . Therefore, we would have to keep T small in 
order to avoid an error that is of the same size as the solution. In spite of 
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this possibility, the scheme is sometimes useful because of its simplic- 
ity. Thus in moving from an ordinary differential equation, we see that 
we must make a scheme that has more consistent accuracy. 

We look now instead for a scheme that is accurate in the equation or 
consistent with the equation to second order in the difference h taking 
k, h of the same order of magnitude. 

We note that a derivative is much more accurately described by a 
difference ratio that straddles equally the point where the approximation 
is made. Thus using 



2k 

for dU/dt the derivative is accurate to second order in k. This can be 
seen as follows. Let W(t) possess a Taylor's series. Then 

W(t + At) = W(t) + W'(t) ■ At + W"(t) ■ (At) 2 /2 + 0((At) 3 ), 
W(t - At) = W(t) - W\t) ■ At + W"(t) ■ (At) 2 /2 + 0((A0 3 ). 

On subtraction, we obtain 

Thus, we are led to a scheme consistent to 0(/z 2 ): 

u j+l -u j - 1 M.-uJ 

_J !_ + a J^±l izl = o 

2k ' 2h 

or 

= U ^ ~ ir^ ~ U !-J (L45) 

Here 'a' is considered as a function of x, t. This algorithm gives us a 
new row of values from neighbours in the two preceding rows. 

Its apparent disadvantage is that it requires the values on two rows 
to start with, but we have only one. The simplest way to get the second 
row is to use the first scheme (ll.44t . This leads to an error of order h 2 in 
the second row. This error is just carried through, but not compounded. 
From the algorithm, we see that the value at the point (i, j) is obtained 
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successively from a pyramid of mesh points and that the algorithm is 
applied j(j - l)/2 times again; but the error is of order h 3 at each step, 
and hence the error is of order j Jr. Since j =* T/h, T being the final 
time, we see that the net error is of order h. 

This is then a consistent scheme but it is still possible for it to be 
unstable, i.e., to have the property that it amplifies the error. 

Stability. (I1.45t (and also (ll.44t ) is a difference scheme with constant 
coefficients. A constant raised to a power plays the role that an exponen- 
tial plays in differential equations. So, we look for solutions of dl.45t 
in the form £'£ J ', i.e., Uj = Substituting this in dl.45b . we obtain a 35 
relation between £ and f . 

^r'-ftf-r 1 ). 

Thus for every £ there are two possible value for £. We are, of course, 
looking for real solutions and we could generate the solutions with real 
But, clearly, if £ and f are complex, then 

(^ ; ' + ?T)/2 and (t't 1 -??)/! V^T. 

will both be real solutions. For j = 0, these solutions are cos(/ arg £), 
sin(/ arg £). If we take \£\ = 1 and put aigg/h = n, n = 0, ±1, 
±2, . . . the corresponding values for j = then become cos(nih) and 
sin(raTz), which can be used to describe any initial data by a Fourier 
series approximation. Let us consider first cos(ra'/z). At an arbitrary row 
j its value, since \g\ - 1, is Re[exp( V— 1 inh)^} where we can take ^ to 
be either root of 

h 

, 2ak i — 

= r V-L • sin(nh). (1.46) 

h 

Note that the absolute value of the product of the roots is 1. So, if one 
of the roots has absolute value different from 1, then there is a root with 
absolute value greater than 1. Suppose j = T/k; then the solution is 
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If the mesh size k shrinks, this solution behaves like exp((Tlog £)/k) 36 
and it also oscillates. The exponential factor goes to infinity as T tends 
to infinity. A small multiple of this solution even, say, of order k , for 
the sake of argument, still goes to infinity. Such a small multiple can 
easily represent an error and hence, the error amplifies as the mesh size 
shrinks, i.e., the scheme is unstable unless \(\ = 1. In that case set 
C = exp( V-T3>), with <I> a real. From ( 11.46b . we then obtain 

/ — / — ak 
2 V- 1 sin <D = -2 V- 1 — sm(nh) 
h 

or 

sin O = sin (to). 

h 

ak 

But this has a solution iff | — | < 1, and hence if this is the case the 

h 

ak 

difference scheme is stable. If | — | > 1, then there is a root £ with 

h 

\£\ > 1 and the difference scheme is unstable. 

What has been established here is the Courant-Friedrichs-Lewy con- 
dition for the particularly simple hyperbolic system of one equation. 

Note. If we look at the characteristics, we see that h/k must be greater 
than characteristic speed for stability. Otherwise, we are trying to eval- 
uate a solution at points whose domains of dependence include points 
from which we are drawing no data. 



Exercise 1.5. Find the stability criterion for <ll.44b . 

We now consider a difference scheme for the nonlinear equatior0 

U t + F x = 0, -oo<x<oo, t>0. (1.47) 

Here U is a scalar and F(U) is smooth. Initial data will be prescribed 
for (OTl) : 

U(x,0) = <D(jc). (1.48) 



We know that in general smooth solutions of dl.47t do not exist for all 
time however smooth O may be; we have to consider weak solutions. 
We recall the definition: 



2 The nonlinear system can be handled inthe same way both formally and numeri- 
cally proviced the speeds of propagation are always distinct. 
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Definition . A locally integrable function U (x, t) is a weak solution of 
Q71 with initial data (Qgt if 

IT (W t U + W x F)dxdt + J W(x,0)®(x)dx = (1.49) 

f>0 

is satisfied for all smooth functions W which vanish for large |x|, t and 
t - 0; we call such functions, test functions. 

We now propose and discuss a difference scheme for getting an ap- 
proximate solution to dl.47t with initial data dl.48t . 

Choose G : M? e — > R, smooth enough, related to F by the reouire- 
ment 

G(u,...,u) = F(u). (1.50) 

21- arguments. 



For k an integer put Ut = U(x + kAx,t) where Ax is step size in x- 
direction; similarly, let At be step size in £-direction. Define 

G(x + -Ax) = G(U- M , U-u2,...Mt) 

and 

G(x - l -Ax) = G(U. e , U. M U^). 
We now consider the following difference analog of (11.47b 

AU AG 

— + — = 1.51 

At Ax 

where 



AU = U(x, t + At)- U(x, t), 
AG = G(x + jAx) - G{x - -Ax). 
It follows from dl.51l > that 

U(x, t + At) = U(x, t) - AAG, 



(1.51)' 
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where A = At/ Ax. 



We claim that the difference scheme dl.51t . as a consequence of 



dl.50t . is 'consistent' with the differential equation (11.47ft in the fol- 
lowing sense: denote by V(x, t) the solution of the difference scheme 
where, we have taken V(x, 0) = ®(x). Here V is defined for non- 
integer multiples t of At, for the sake of convenience, as equal to V(x, f), 
t' = [t/At]At. Ofcourse, V depends on Ax, At. 
Then the following theorem holds. 

Theorem. Assume that as Ax, At — » 0, V(x, t) — * U(x, t) boundedly a.e. 
Then U(x, t) is a weak solution of \1.47\ with initial data M.48\ . 



Proof. Multiply dl.51t throughout by At and by any test function W and 
then integrate with respect to x to obtain 



r av r ag 

W(x,t)—dxAt+ W(x,t)— 
J At J Ax 



dx At = 0. 



Now sum over all t which are integral multiples of At and carry out 
summation by parts in the first integral; we obtain 

v-i f w(x, t - At) - w{x, t) r 

2^ I — — ^ K -^V{x,t)dx At - I W{x, 0)®(x)dx+ 

r AG 
2^ J W{x,t)—dx At = 0. 



t>0 

+ 



In the last integral replace x by x - - Ax in first term and by x + - Ax in 
second term; we finally obtain 



2/ 



W(x,t-At)-W(x,t) 

V(x, t)dx At- 



At 

t>0 



-J W(x, 0)®(x)dx 

^ r W{x + 1/2 Ax) - W(x - 1/2 Ax) 

- > GdxAt, 

V J Ax 
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where G stands for G(V\, . . . , V%i), V\,..., Vit denoting values of V at 
2( points which are distributed symmetrically around (x, t) and have dis- 
tance Ax from each other. If V — > U boundedly a.e. as Ax, At — » so 
do Vt , . . . , Vze an d therefore 



G(V U . . . , V 2 {) -» G(f/, ...,£/) = F(U) by (L5Q) . 
The proof is complete. □ 

The real difficulty is to find when V(x, t) U (x, t) boundedly. 

We turn to the problem of choosing G and minimizing the truncation 40 
error. Let U(x, t) be an exact smooth (C 2 is enough) solution of dl.47t . 
It will then satisfy difference equation (1.51)'| only approximately; the 



deviation of right side from the left side of |(1.51) / ] is called truncation 
error. It is easily seen that, in view of (ll.50l l. then truncation error is 
0(A? 2 ). We shall now show, by taking i — 1, that G can be so chosen 
that the truncation error is 0(Ar ). Let 

U(x, t + At) = U(x, t) + At U, + ^(At) 2 U t! + 0(At 3 ) (1.52) 

be a Taylor series up to terms of second order. 
From (11.471) . we obtain the following: 

U, + F X = 

2 (1-53) 
U tt + (A 2 V x ) x = 0. 

The second of these equations follows from the calculation with 
dF/dU = A, 

U„ = -F xt = -F tx = -(AU t ) x = -{AF X ) X = -(A 2 U X ) X . 

What is significant is that all t derivatives are exact x derivatives and 
therefore can be approximated by exact x differences. Substituting 
(ll.53t in (11.521) . we obtain 

U(x, t + At) - U(x, t) + (AtF + j(At) 2 A 2 U x ) x + 0(At 3 ). (1.52)' 
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Comparing (1.51)' and (1.52)'[ we see that the truncation error is 0(A? 3 ), 
iff 



AG 



1 



— = (F + -At A z U x ) x + 0(At z ). 
Ax 2 

41 From this, we can easily determine the form that G must take. 

Theorem. The truncation error in the difference scheme (1.51)' is 0(A? 3 ) 
iff 

F(a\ + F(h\ 1 

(1.54) 



G(a,b) = F(a)+ 2 Fib) + 1 -AA 2 .(b-a) 



plus terms which are 0(|a - b\ 2 )for (a — b) small. 



The quantity A 1 in (ll.54t shall be taken as l/2{A 2 (a) + A 2 (b)} for 
the sake of symmetry more than anything else; any other choice would 
make a difference that is quadratic in (a - b). 

Denote the function in dl.54t by G ; any permissible G can then be 
written in the form 



G = G + ^Q(a,b)-(b-a) 



(1.55) 



where Q(a, b) vanishes for a = b. Substituting ( 11.551 ) in (1.51)'[ we see 
that 

U(x, t + At) = U{x, t) + AA'F + ^A 2 AA 2 AU + -AAQAU, (1.56) 

where A' - l/2[7\Ax) - T(-Ax)] and A = T(-Ax) - T{--Ax), T(s) 
being the shift operator of the independent variable by an amount s. We 
shall call Q the artificial viscosity. 

The difference equation (11.561) expresses the value of U at time t+ At 
as a nonlinear function of U at time t; we shall write this as 



U(t + At) - NU(t). 



(1.56)' 
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The value of the solution of the difference equation at some later time 
kAt is obtained from the initial values by the application fc-times power 
of the operator ^V. 
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Our aim is to show that the difference scheme dl.56t is convergent 
if the size of A is suitably restricted. In the case of linear equations, it is 
well-known and easy to show that convergence is equivalent to stability 
defined as the uniform boundedness of all powers N k of N with some 
fixed range kAt < T, a problem which we have studied previously. In the 
nonlinear case, following Von Neumann the convergence of the scheme 
would depend on the stability of the first variation of the operator N. 
The first variation of N, by definition, is a linear difference operator 
with variable coefficients; Von Neumann has conjectured that such an 
operator is stable iff all the localized operators associated with it, i.e., the 
operators obtained by replacing the variable coefficients by their value 
at some given point, are stable. 

We content ourselves with: 

Theorem. If the Courant-Friedrichs-Lewy condition 



is satisfied, the difference equation M.56\ satisfies Von Newmann 's con- 
dition of stability that the linearized equation is stable. 

Proof. The first variation of the operator ,/V can be easily computed and 
is given by 



where A', A are as before, and 0(Ax) denotes an operator bounded in 
norm by |Ax|, provided we are perturbing in the neighbourhood of a 
smoothly varying solution, i.e., one where neighbouring values differ 
by Q{Ax). In this case, the influence of the additional viscosity term is 
§{Ax); see Remark ffl below. 

To localise the operator dl.58t . we replace A by its value at some 
point. After making a Fourier transform, the operator A' becomes mul- 
tiplication by i sin a, and the operator A multiplication by 2i sin( 1 /2)a, 
so that ( 1 /2)A 2 becomes multiplication by cos a - 1 ; here a = £Ax, £ 
being dual variable. Hence the amplification function of the operator 




Imax 



(1.57) 




(1.58) 



T38lis 



/ + iAA sin a + A 2 A 2 (cos a - 1) + 0(Ax). 



(1.59) 
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1. The Traffic problem and a first order nonlinear equation 



Since the eigenvalues k of A are real, the eigenvalues v of the matrix 
dl.59l > are given by 

|v| 2 = (1 - k 2 {l - cos a)) 2 + k 2 sin 2 a + O(Ajc) 

= 1 - 2k 2 {l - cos or) + k 4 (l - cos cr) 2 + k 2 {\ - cos 2 a) + 0(Ax) 
- 1 - (k 2 - k 4 ){\ - cos a) 2 + 0(Ax). 

By our assumption dl.57t \k\ < 1. Thus |v| < 1 + 0(Ax). 

The proof of completed. □ 

Remark 1. We have seen above that the quadratic terms, Q(a, b), ap- 
pearing in a representation of G{a, b) influence neither the order of the 
truncation error nor the stability of the scheme at points where the so- 
44 lution varies smoothly. The terms do influence, however, at the points 
where solution varies rapidly, e.g., across a shock. For a detailed analy- 
sis, we refer to Lax and Wendroff [24]. 

As a second example, we consider the Lax-Friedrichs scheme for 
(11.47b . We add a dissipative term eU xx with e - 2Ax 2 /At. This scheme 
is given by 

U'] +l = U'j - J±L[F(U] +1 )- F{U n hl )] + ±(U n j+1 + U n j - 2Up, (1.60) 

where U" abbrevates U(jAx,nAt). We want to establish convergence 
via the contraction mapping principle. We write dl.60t as 

Then T maps sequence jf/;} ._ to a sequence {r(£/)} . according 
to 

{T(U)}j = Uj - ^[F(Uj +1 ) - F(U hl )] + l -(U j+l + U H - 2Uj), 

Let l x be the space of all summable sequences {£/,} _ with usual or- 
dering: if U, V € € l say U < V if Uj < Vj for all j. Let a < < b and 
put 

C - {U € e 1 : a < Uj < b for all j}. 



1.7. Computing methods 
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Assuming F(Uj) - F(U-j) — » as y —> oo it can be seen easily that 

DO OO 

y=-oo ;=-oo 

Thus T is integral preserving. It is also easy to see that if At / Ax\F' (Q.)\ < 45 
1, a < D < b, then T is order preserving, i.e., U < V. 

T(U) < T(V). Then it follows, by the following theorem, that T is 
also a contraction. Thus the scheme is convergent. 

Theorem . Let Q. be a measurable space with a positive measure and 
T : L\Sl) -> L\Q) satisfy 

/" /' 

Let C L'(Q) be such that whenever f, geC, max(/, g) € C. Then the 
following statement are equivalent. 

i) /, g £ C, f < g a.e. T(f) < T(g): order preserving. 

ii) / (Tf - Tg) + <J(f- g) + , fygeC where r + = max(r, 0). 
n n 

iii) J \Tf -Tg\< J \f - g\, f,ge C: contraction, 
a a 

For proof, we refer to Crandall and Tartar 



Chapter 2 



One Dimensional gas 
dynamics 



2.1 Equations of motion 

These equations of motion completely characterise smooth movement 46 
of a fluid. They express the physical laws: 

i) Conservation of mass, 

ii) Conservation of momentum and 

iii) Conservation of energy. 

We first consider conservation of mass. For this, the rate of change 
of mass in any volume element V of the fluid is balanced by the flow 
across dV, the boundary of V. If n denotes the outward unit normal to 
d V, then the normal component of velocity across dV is n ■ u, where u = 
(wi,W2,W3) is the velocity vector of the fluid (we are mainly interested 
in 1, 2 or 3 dimensional flows). Then the net flow across the boundary 
between two times t\ , ?2 is 




h dV 
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2. One Dimensional gas dynamics 



where p is the density of the fluid and dcr denotes the surface measure. 
This must be balanced by the change in total mass between the times t\ , 
t2- Hence 



[ J pdV]'* = -J j pin- u)do- dt, 



h dV 

or taking the limit as t\ — > t%, 



^- J pdV + J p(n ■ u)do- = 0. 



V dV 

47 Using divergence theorem this can be written as 



j t J pdV + J div(pu)dV = 0. (2.1) 

v v 



In a similar way, the equation for the net change in the i component of 
the momentum is 

d f C 

— I piijdV + I [pui(n ■ u) + pni\dcr = 0. (2.2) 

v dV 

The first term is the rate of change of total momentum inside V, the 
second term is the transport of the momentum across the boundary and 
the third is the rate of change of momentum produced by the pressure p. 
Here, we are neglecting other forces such as gravity, viscosity, etc. 

The total energy density per unit volume consists of the kinetic en- 
ergy p\u\ 2 /2 of the motion of particles plus the internal energy pe of the 
molecular motion. For energy balance, we then have after neglecting 
heat conduction, viscosity, etc. 

j t J (p|«| 2 /2 +pe)dV + J~{(p\uji 2 + pe)n • u + pn- u\dcr - 0. (2.3) 

v av 

The first term in the surface integral is again the contribution from en- 
ergy transport across the boundary and the second term is the rate of 
work by the pressure p at the boundary. 



2.2. Thermodynamical relations. Entropy: 
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If discontinuities are allowed in the flow, we have to work with d2.lt 
- (12.31) (In fact they are not quite enough). But if all the quantities are 
smooth, we can different under the integral sign. We then obtain the 
differential equations 

+ div(p M ) = 0, (2.1(a)) 
ot ~ 

d d 

—(pud + div(pwjM) + — = 0, i = 1, 2, 3, (2.2(a)) 

Ot OXi 

—(p\u\ 2 /2 + pe) + div[pu(p\u\ 2 /2 + e + -)] = 0. (2.3(a)) 
ot ~ ~ ~ p 
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2.2 Thermodynamical relations. Entropy: 

Consider the quantity de + pd(l/p). When a small energy is added to 
a mass of gas some of the energy is the work of changing the volume 
l/p to l/p + d(l/p). This energy is pd(l/p); the rest de is heat put into 
the system. Since the quantity de is a perfect differential, there exist 
functions S(p,p), T(p,p) such that 

de + pd(l/p) = TdS (2.4) 

T is the absolute temperature and 5 is the entropy. The relation (12.4ft 
may be viewed as a way of defining temperature and entropy up to an 
arbitrary function. 

Suppose we treat p, r = p~ l as the independent variables, then we 
will have 

de = -pdr + adp + (3dr, 

where a, ft are given functions and we want to express adp+fidr as TdS , 
where T, S are functions of p, r. We want to show this representation 
is essentially unique. First from the compatibility relations for de, we 
have 

(rP+fi) P = oiT, (*) 
and from the one for S we have 49 
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{a/T) T = (fi/T) p 

or 

a(-) T -/3(-) p + (a T -p p )-=0 

or 

a(-) T -J3(-) p - - = 0, from(*). 

This is a first order equation for \/T which we may solve provided a 2 + 

1 ± 0, which is a reasonable thermodynamic assumption. Once T is 

determined S is found from dS = (a/T)dp + (J3/T)dT. With respect to 

uniqueness, suppose we have two such representations, T, S and T',S'. 

Then TdS - T'dS' or S T /S P - S' T /S' p and hence F{S,S') - for 

some F or with an assumption of monotonicity S' - s(S) and then 
T 

T = Thus the temperature and entropy are uniquely determined 

ds/dS 

up to a scaling factor ^(5). 

Suppose, we know for some medium that e is a function of t and 5 , 
where t - lip, specific volume. Then we see that p is also a function 
of p, S . We assume 

p = f(p,S) or p = g(T,S). 

It is a fundamental property of almost all media that, entropy remaining 
constant, pressure is an increasing function of p or equivalently decreas- 
ing function of t. Thus 

f p > and g T < 0. 

50 For any value of S, the function g(r,S) is generally convex w.r.t. t. 
Henceforth, we shall assume this: 

g TT (r,S)>0. 



Alternative forms of equations of motion: 

Introducing the operator 

D_ _ d_ d_ 

Dt ~ di + U] "d^j' 



2.3. One dimensional flow 
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the equations ( |2. l(a)| > - ( |2.3(a)| i can be written as 




(2.5) 



(2.6) 



(2.7) 



or using (12.51) . (12.7ft can be written as 




De p Dp 
~Dt~~fi~Dt 



Then from the thermodynamic relation ( 12.4ft . we obtain 



T = 0. 

Dt 



(2.8) 



This is to say, entropy remains constant following a particle. Flows 
satisfying ( 12.8ft are called adiabatic. It follows from that, if the fluid 
initially has uniform entropy then entropy remains constant throughout 
as long as the flow is continuous. In such a case p - f{p). Such flows 
are called isentropic. 

2.3 One dimensional flow 

The equations reduce to 



The second equation in ( 12.91 ) can be written, using the first equation, as 51 




(2.9) 



pu t + puu x + p x = 0. 
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The third equation follows from d2.8t . It is convenient to use p, u. Then 
the system i2.9l may be written as 

p t + up x + pc 2 u x - 0, 

pu t + p x + puu x = 0, (2.9(a)) 
S t + uS x = 0, 

where c 2 = (dp/dp)s = constant; c is called the sound speed. System 
( |2.9(a)| ) is a typical nonlinear hyperbolic system. 

Another such ta system, the Lundquist equations, occurs in magneto- 
hydrodynamics. It involves coupling Maxwell's equations with the clas- 
sical equations of gas dynamics. 

Let E, B stand for electric and magnetic field vectors and let J, u 
denote the electric current and the flow velocity. If p is pressure and 
c 2 = dp I dp is the sound velocity, then the equations are 

dB 

— + (V-u)B-(B- V)h - 0, 

dt 

du 

P^r + Vp + B x (V x B) = 0, 

dt 

c 2 lp ■ ^ + c 2 V • u = 0, 
dt 

dt 

52 Another related hyperbolic system occurs in the theory of shallow water 
where (I2.5t and (I2.6t hold with the depth playing the role of density p 
and pip 2 is constant. 

We turn to the study of the system ( |2.9(a)| i. We consider isentropic 
flow so that the third equation drops out. Adding and substracting the 
other two equations, we arrive at 

p t + (u + c)p x + pc{u t + (u + c)u x } = 0, 
p t + (u - c)p x - pc{u t + (u - c)u x ) = 0. 

We choose the characteristics C+, C_ so that the derivatives are along 
them. Thus for 

r dx 

C + : — - u + c, 
dt 
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dx 

C_ : — = u — c 
dt 

and then 

dp du 

— + pc— = on C+ 

dt dt 

dp du 

— — pc— =0 on C . 

dt dt 

Integrating these two equations, we obtain 

€{p) + u = constant on C+, 
£{p) - u = constant on C_, 

where 

up)- r cip)dp 



p 

The quantities €{p) ± u are called Riemann invariants. 

Exercise 2.1. Reduce the system ( |2.9(a)| i (omitting the last equation) 53 
for slow speeds nearly constant density to the wave equation (in u or in 
p). Show also that discontinuities in derivatives propagate with sound 
speed. 

Exercise 2.2. Find the speeds of propagation of the one dimensional 
Lundquist equations. 



Simple Waves. Piston Problem: In addition to being isentropic, if 
the flow has one Riemann invariant constant throughout, the solution is 
called a simple wave. Such a wave lies next to a constant state since 
all the characteristics of one kind issuing from the constant state carry a 
constant Riemann invariant. 

To illustrate, how such waves can be produced, we consider the pis- 
ton problem as a basic model. Consider the waves produced by the 
movement of a piston at the end of a tube and gas at rest in a constant 
state ahead of it. Provided shocks do not appear, the waves carry a con- 
stant Riemann invariant from the constant state ahead. 
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Initially, the gas has velocity u = 0, sound speed c Q and S - S in 
s > 0, t = 0. The piston path is itself a particle path and the flow is 
isentropic. 

Since u - c < u the C characteristics start on the x - axis in the 
uniform region (see figure 2.1). On each C_ 

£{p) — u — constant; (2. 10) 

54 Since u = initially, we conclude that the constant is £{po), where po 
corresponds to the initial density. Since this constant is the same for C_ 
characteristics, we conclude that £{p) - u is the Riemann invariant that 
is constant throughout. We now use the other characteristic to find the 
other quantities. 




u = 0, c = c , S = S 

Fig. 2.1. 



For those C + which originate from the x-axis, we see that 

€{p) + u = Ap ) (2-11) 

From (I2.10I ). ( 12.1 U . we conclude that u = 0, c = Co in the region covered 
by these C + characteristics. Let C% separate the C+ characteristics that 
meet the piston from those which meet the x - axis. 

Since the image of the flow in (p, u) space lies entirely on the curve 
i'{p) - u = constant each point on the curve represents a curve in the 
flow. At that point the value of tip) + u determines u,p and hence the 
slope dx/dt - u + c, of C+ which is constant along that characteristic. 
Hence 



2.4. Shock conditions 
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p, u are constant on each C + : — - u + c. 

at 

To complete the solution, we should use the boundary conditions, 55 
given on the piston. Let the piston path be given by x = X(t). Then the 
boundary condition is 

u - X(t) on x- X(t). 

With this, we can obtain the complete solution: 

x = (X(t) + c) (t - t) 
u = X(t) 

Since the C + characteristics are straight lines with the slope dx/dt 
increasing with u, it is clear that the characteristics will overlap if u in- 
creases on the piston i.e., if X(t) > for any t (this is the case when the 
piston accelerates into the gas). This is typical of nonlinear breaking, 
(see p8) and shocks will be formed. We need to reexamine the argu- 
ments leading to the constancy of entropy and of one of the Riemann 
invariants. But for motions with X(t) < 0, we have constructed a com- 
plete solution. 

2.4 Shock conditions 

We suppose the flow is one dimensional and allow jump discontinuities, 
denoted by [■], in the flow. We must work with the equations d2.lt - 
(I2.3I) . Let the region V collapse on a slit around the discontinuity and 
we find the discontinuity conditions: 

-U\p] + \pu] = 
-U[pu] + [pu 2 + p] = 
-U[p 2 /2 + pe] + [(pu 2 /2 +pe + p)u] = 0, 

where U is the discontinuity velocity. We denote by subscripts (0), (1) 56 
the two states. As a particle moves across the discontinuity, it moves 
from the front of the discontinuity across the discontinuity to behind the 
discontinuity. The entropy must increase in this direction. 
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For convenience, we introduce the relative velocities 

Vi = U - Uj, i = 0,1. 

In the velocity frame where the discontinuity is at rest, since the equa- 
tions are independent of the coordinate system, we obtain from the 
above 

PqVq = pm - m, (2.12) 
where m is the mass flux through the surface; 

PoVq + Po = P\v\ +P1-P, (2.13) 

where p is total momentum flux and finally the energy flux condition: 

m(vl/2 + e + pqTq) = m(y\l2 + e { + p\Ty). (2. 14) 

According to the second law of thermodynamics entropy can only in- 
crease. Hence across a discontinuity 

S <Si 

or 

mSo < mS i (2.15) 

Two two types of discontinuity surfaces are distinguished by the cases 
m = and m ± 0. The case m # corresponds to mass flowing across 
57 and the discontinuity surfaces are called check fronts; the case m = 
corresponds to a contact or slip surface. 

We first consider the case of a shock, m + 0. Then from equation 
(12. 14b . we obtain 

1 2 1 2 

-v + e + p T = -Vj +ei + pm. 
If we introduce enthalpy, i = e + pr, we then obtain 

vl/2 + i = v\/2 + h (2.16) 
If we use (12. 1 2b and <I2. 13b . we obtain 



ToiPO- P\) = Vq(V\ - V ), 



2.4. Shock conditions 
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TlOo - Pi) = Vl(Vl - V ) 

and hence adding these two equations, we obtain 



(to + T\) (pi - p ) = Vq - v\ (2.17) 



Equation (I2.16t becomes 

Oi - Po) (to + ti)/2 = h - k, (2.18) 

or, since i — e + pr, 

(j -T 1 )(p l +po)/2 = e l -e . (2.19) 

Since d2.18t and ( I2.19t involve only thermodynamical quantities, they 
are particularly useful. They were first studied by Hugoniot and d2.19t 
is called the Hugoniot relation. The relation ( I2.19t can be intepreted by 
stating that the increase in energy across a shock is due to the work done 
by the mean of the pressures in performing the compression. 

The most common example considered is a poly tropic gas, with p - 58 
A(S)p y . We note the following relations: (refer Serri 



ropr 



log A = S/c v , e = pliy - l)p, c 2 = yp/p, 

where c v is the specific heat at constant volume. It is often useful to use 
formulas for polytropic gases with the parameter 

M = ^^ 

which is the Mach number of the shock relative to the flow ahead and is 
a useful measure of strength. Using ( 12.121 1 and ( 12.131 1 and the definition 
of M, we arrive at the following: 

■ mi - wo . = 2(M 2 - 1) 
c ( r +l)M' 
Pi/po = (y+l)M 2 /{(y-l)M 2 + 2}, 



'See, Fliigge, S. (Ed.), Encyclopedia of Physics, VIII/1, Springer- Verlag (1959). 
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ip\ ~ P0)/P0 = 2y(M 2 - l)/(y + 1), 

{2yM 2 - (y - 1)} 1/2 • {(y - 1)M 2 + 2} 1/2 

c i/ c = ; 7777 • 

(y + 1)M 

Introducing an alternative strength parameter 

Pi - Po 
z = , 

Po 

the above set of equations can be written as 

M -_{ l + ^.) m , (2 . 20) 

(2.21) 



2y 

U\ — Uq 



Pl/P0 = { 1 + <I±S} / { 1 + <I_i)£}, ( , 22) 

(d + z)(i + ti£)|" 2 
"*•"{ 1 + ^ ) <2 " 23) 

59 From the relation for 5 , obtained from A, and the relation (12.22b . we 
obtain 

f(l + z)(l + ^ )z ) 
(5 1 - S )/c v - log 7 [ (2.24) 

Since 5 1 > 5o across a shock, we obtain, as is easily seen by expanding 
that z > 0. Hence p\ > po, and from the above relations, we then obtain 

Pl > po, C\ > Cq, U\ > Uq, M > 1. 

From ( I2.20I I. it is clear that U > uq+cq. From the relations (12.20b . (I2.21t 
and J2.23l >. it then follows that, 

u\ + c\ > U. 

Thus a poly tropic shock is always compressive with p\ > po and it is 
supersonic viewed from ahead and subsonic from behind. 

These facts are true quite generally, but we will discuss this later, 
when we study the Riemann problem. 



2.5. Contact Discontinuities 
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2.5 Contact Discontinuities 

If the mass flux m through the surface of discontinuity is zero, then 
Vi = vo = 0, so that ui - uq - U. Then we have po = p\ from 02.13b 60 
and (12.141 1 is automatically satisfied. Such a discontinuity surface as 
indicated before, is called a contact or slip surface. 

The flow velocity is continuous across the contact surfaces in one 
dimensional flow, but in higher dimensions the tangential component of 
the velocity vector may suffer a discontinuity across a contact surface, 
while the normal component relative to the surface is always zero. For 
details see, e.g., Courant-Friedrichs [6]. 

2.6 Shock Reflection 

A simple example of determining a flow with shocks is provided by 
the reflection of a shock from a wall which can also be solved exactly. 
A shock with a state behind of prescribed velocity hits a wall and is 
reflected. We seek the pressure after reflection. 

Let the subscripts 0, 1 refer to the states ahead of and behind the 
incident shock, and subscript 2 refer to the state behind the reflected 
shock, see figure 2.2. 
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Fig. 2.2. 



61 If the shock strength of the incident shock is zi = ip\ - Po)/Po, the 

state (1) can be determined by the relations (12.201) - <!2.24t . If 

Zr = (P2 ~ Pi) I Pi 

is the strength of the reflected shock, we obtain, with suitable change in 
the sign of the velocities since the reflected shock travels in the opposite 
direction to the incoming shock, from (12.211) that 

,"l-M2, Zr 



Next to the wall, the gas must be at rest and hence U2 = «o = 0. Now, u\ 
and ci can be found in terms of zi- After doing this, we obtain 

Zi Zr 



y{(i+z/)(i + ^)} 1/2 rU + ^) 1/2 

This is a quadratic in zr and it is easily seen that the relevant solution is 

Zi 



Zr 



, , (y-i)z; 

1 + 2y 



2.7. Hugoniot Curve. Shock Determinacy 
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For weak shocks zi — > and hence z« - z/ and 

p 2 - po - 2(pi -po). 

So, for the acoustic case, the pressure is doubled as is well known. Also 62 
for strong shocks zi — > 00 ■ This implies z# = and therefore 

(y-l) 

P2 3y- 1 _ 
Pi 7-1 

for y = 1 .4. Hence there is a large gain in pressure after reflection. This 
phenomenon is even more striking in spherically symmetric flows when 
a shock is reflected at the center of symmetry. 



2.7 Hugoniot Curve. Shock Determinacy 

We recall the Hugoniot relation ( 12.191 : 

e\ - e = (t - Ti) (pi + po)/2. 

We regard all quantities such as energy e, entropy 5 etc. as functions of 
t, p. We define the Hugoniot function 

H(t, p) = e(r, p) - e(r - po) + (t - t ) (p + po)/2. 

If (to, po) is fixed, the graph of the points (r, p) which satisfy H(t, p) = 
is the Hugoniot curve. As we shall see with p > po it represents all 
possible states that can be reached with (tq,po) ahead of the shock. For 
p < po the curve represents the states that can be ahead of (to, po)- The 
big advantage in considering this curve is that it involves no velocities. 
For polytropic gases, we have 

_ PT _ PT{\ - p. 2 ) 

(y-l) 2p 2 ' 



where 
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Hence 



2)/H(t, p) = ( T - h 2 tq)p - (t - h 2 t)pq. 



Hence the Hugoniot curve, specifically the Hugoniot curve with center 
(to > Po)-> is a rectangular hyperbola with the left asymptote 



t = fX 2 T = T min > 0. 



(See figure 2.3) 




, (to - t)(p - po) = ("l - w ) 2 
Front state with state (0) behind 

(T - T )(po ~P) = («1 - M()) 2 

Bock state with state (0) in front 



T 



T 



Fig. 2.3. The Hugoniot curve is the heavy curve. 

We shall see that under wide conditions on the Hugoniot curve the 
state (r,p) that can lie ahead or behind (ro,po) an d connected by a 
shock, can be found exactly. 

The precise statements are: 

1. The state (0), i.e., in which (to, po) is given, and the shock velocity 
U determine the complete state (1) on the other side of the shock 
front. 

2. The state (0) and the velocity u\ determine the speed of the shock 
front and the complete state (1) if it is specified whether the state 
(0) should be ahead of or behind the shock front. 



2.7. Hugoniot Curve. Shock Determinacy 



57 



3. The state (0) and the pressure p\ determine the speed of the shock 
front and the complete state (1). 

To prove these statements, we make the following assumptions on 
the Hugoniot curve H(t, p) = with center (tq, Po)'- 

(HI) Along the Hugoniot curve, the pressure varies from zero to infin- 
ity and the value of r exceeds r m j n . 

dp 

(H2) The Hugoniot curve is strictly decreasing, i.e. — < along the 

dr 

curve. 

(H3) Every ray through (to, Pq) intersects the Hugoniot curve at exactly 

d 2 p 

one point and at (tq,Po), — 7 > 0. 

dr 1 

(All three conditions are satisfied by the polytropic gases). 

We are now in a position to prove statements (1) - (3) made above 
provided we assume pressure increases across a shock. It follows from 
the shock condition (I2.12t and (12.131) that 

2 2 2 Pi ~ P0 ,r. OC x 

-m = p v (2.25) 

So, to find (T\,pi), we need to find the intersection of the curve with 65 
the line through (To,p„) and with slope —nr. Then (H3) assumes there 
is just one such intersection. The velocity u\ can then be found through 
\2.Yl\ . U] = U — niT\. This proves statement (1). 

As far as second statement is concerned, it can be easily derived 
from (12. 1 2ft and (12. 1 3I > that 

~(tq - Ti) Oo - Pi) - (vi - v ) 2 . (2.26) 

From the data given, {u\ - uq) 2 - (y\ - vq) 2 , is known. Hence to find 
(ti,/?i) it suffices to find the intersection of the hyperbola 

(to - t) (po ~p) = -(vi - v ) 2 

with the Hugoniot curve. The slope of the hyperbola is m 2 > 0. Hence 
from (H2), it follows that there are just two such intersections, corre- 
sponding to the two possibilities that the state (0) lies ahead of or behind 
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the shock front (see figure 2.3). The shock velocity and the state (1) can 
then be found easily using shock conditions. 

To prove statement (3), assume the state (0) and p\ are given. The 
assumptions (HI) and (H2) assert that there is exactly one t\ such that 
H(t\ ,pi) = 0. The other quantities to determine the state (1) completely 
are then found essentially as above. 

We shall discuss a few more qualitative statements about the shock 
transition using the Hugoniot relation; in particular, some of the proper- 
66 ties of the shock transistion were already discussed in an earlier section 
for polytropic gases. The main thing we are going to prove here is the 
following: 

The increase of entropy across a shock is 
of the third order in the shock strength 
and the shock is compressive. 

By the shock strength, here, we mean one of the quantitites 



We can consider, because of (H2), the Hugoniot curve as p = G(t); 
in particular, we can consider t as independent variable. Now along the 
Hugoniot curve dH = 0. So, 



Pi-Po, Pi- Pa, or |vi — vol- 



2de + (t - To)dp + (p - po)dT - 0. 



But de + pdr - TdS and therefore we obtain 



2TdS - (p - po)dT + (t - To)dp = 



(2.27) 



and hence 



dS=0 at (to,po)> the center. 



(2.28) 



Differentiating H2.27II . we obtain 



2d(TdS) + (T-T )d 2 p - 0, 



(2.29) 



and hence again at the center (tq, Po) 



d(TdS) = dTdS +Td 2 S -0 



2.7. Hugoniot Curve. Shock Determinacy 
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and therefore from d2.28t 

d 2 S = 0. 

Thus the change in entropy is at least of third order. We show it is of 67 
third order. Differentiating A2.29I ). we obtain, at the center {tq,pq), 

2Td 3 S + drd 2 p = 0. 

Thus since d 2 p/dT 2 > at (jo,po), we have 

d 3 S 3; if dr $ at (r ,po) 

which proves entropy change is of third order exactly. 

Furthermore, since the entropy must increase, we must have dr < 
which means the shock is compressive and hence the upper branch of the 
Hugoniot curve represents states behind the state (to, po) as was asserted 
earlier. 

Exercise. Piston at uniform speed: It is a simple exercise to find the 
flow if a piston is moved with uniform speed into a gas at rest. The 
speed of the flow behind the shock is that of the piston and so we are in 
case 2. 

Remark. We need the notion of vorticity in three dimensions. The vor- 
ticity is defined by 

oj = curl u. 

Claim. The change in vorticity across a shock is also of third order. In 
a steady flow the vorticity vector oj and velocity vector u satisfy 

2 1 
V(|m| /2) + {u> x u) h — Vp = constant. 

~ ~ P 

Note that if co - 0, i.e., if the flow is irrotational, we have Bernoulli's 
law. Taking a scalar product with a vector \£\ = 1, lying in the shock 68 
surface and letting [•] to denote difference across the shock, we have 

[2t(| H |2 /2 + pip) - p£(l/p)] - [f • (u x co)] = 0, 
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where the differentiation w.r.t. arc length along the shock is denoted by 
d/ds; but since 



and 

we obtain 



\u\ 2 /2 + p/p] = -[e] 



de + pd(l/p) = TdS, 



dS 

[T-] = [<o-(Zxu)l 
ds - 



Thus the change in oj is of third order in the shock strength if oj = on 
the side since we may choose co-ordinates so that u is locally normal. 



2.8 Riemann Problem 

We now turn to another important initial value problem, the Riemann 
problem. It is also referred as the shock tube problem. It is important 
both theoretically as we shall see and because of its practicability; it is 
the main device for producing fast chemical reaction fronts. The Rie- 
mann problem can be stated as follows: 

Given two states they can always be connected 
by a "fan wave" consisting of a centered 
rare-faction wave, a shock and a discontinuity. 

The two different states will be separated by a thin diaphragm up to 
69 time t = 0. Then the diaphragm is instantaneously removed and we have 
to find the flow. Without loss of generality we can always assume that 
ur = and pr < p^. The subscripts R, L refer to right and left states 
respectively. 

We first define a state behind or SB curve associated with a state 
(tr,Pr). It consists of two branches, the upper branch for r < tr is a 
Hugoniot curve and for t > tr it is a curve p = p(r) at constant entropy 
and corresponds to a centered rarefaction wave leading from (tr, pr) 
and on which u + £{p), see section 1231 is a constant. 

Any point on the SB curve, thus defines a unique transition from 
(tr,Pr) to a new state by means of a shock or a rarefaction. For the 
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right side, these must be facing so that particle paths enter them, i.e., to 
the right. 

Similarly, we define a state behind or SB curve for (tl,Pl) in the 
same way again remembering to face the appropriate waves, i.e., to the 
left. 

We also define a map from the 5 B curve of (tr, pr) to that of (jl, pi) 
on the lines p = constant. Since the SB curves are easily seen to be 
monotonic, this map is invertible and is represented by t(t) if r lies 
on the SB curve of (tr,pr) and by f _1 (r) if r lies on the SB curve of 
(r L , pl). See figure 2.4. 

The solution of the Riemann problem is found in the (r, p) plane by 
connecting (j R , p R ) to (t l , Pl) by at 




\ (tl,Pl) 
— \ 



*~ T 



Fig. 2.4. 



most two in between states, each representing states behind (tr, pr) and 
(tl,Pl) respectively and connected by a slip or contact discontinuity 
which appears as the horizontal line p = constant. The velocity u in the 
two states must be the same. From the right hand side, it is determined 
either by the shock condition 
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[Or -p)(j- tr)] 



1/2 
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using ur - or by the rarefaction formula 



u 



= -Mtr) + Mr) 



and on the left side by 



u-u L = -[{p L - p){t- t l )] 



or 



u-u L = £ l (tl) ~ trfj)- 



We note that these formulas are continuous at t = tr and r = tl re- 
spectively. We further note that we connot "add in" rarefaction waves 
or shock waves as part of the states to be connected because we cannot 
match the velocities. 

To solve the Riemann problem we have only to show that we can 
always find a horizontal segment where the values of u on the intersec- 
tions with the two S B curves are the same. 

But for p — > oo these intersections correspond to the shock portions 
of the S B curves, the two values of t are approaching their minima and 
the difference between u on the left S B curve (ul) and on the right (ur) 
satisfies 



On the other hand as t — > oo a line p - constant intersects the two 
rarefaction sections of the S B curves. Both Ir(t) and ^l(t) — > since 



ul - uR — > -oo. 



p 




o 



Hence 



u L - ur -> u L + £ l {tl) + Mtr)- 



Hence if 



u L + 4(tl) + Mtr) > 



there is always a solution. On the other hand if 



u L < -Mtl) ~ Mtr) 
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we have two completed rarefaction waves and there is a vacuum in be- 
tween. 

Thus every Riemann problem can be solved. The configurations are 
as follows: 

If ul lies in the interval 

I-Mtl) + Mf(T R )), {(PR - p(t-\t l )) (f-\r L ) - r R )} 1/2 ] 

we have a rarefaction from the left and a shock from the right. This 
includes the case ul < since t(t«) < Tj,. If ui lies in 

[-h(T L ) ~ MTR), 4(f(T R )) - 4(T L )] 

we have two rarefactions possibly with a vacuum and if 

UL > {(PR ~ P(f-\T L ))) (f-\T L ) ~ T R )} 1 ' 2 

there are two shocks. 

It is also quite easy to show that d(ui - u R ) > provided the Hugo- 
niot curves are star-shaped, which shows that within the class of fan 
waves our solution is unique. 

However, full uniqueness follows only by using a contraction theo- 
rem, see, for example, Oleinik [34] or Keyfitz [20]. 



2.9 Solution of initial value problem 

It was originally proposed by Godunov that the initial value problem 
(IVP) should be solved approximately by considering the initial data 
as approximately piecewise constant and then solving a set of Riemann 
problems, each by what we shall call a fan solution. 
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A t 




a; 



Fig. 2.5. 



This solution is considered up to time At where At /Ax < 1/25 
where S is the maximum of shock speeds or characteristic speeds that 
occur. This cuts out intersections. At the next time t we have to replace 
the resulting data again by a piecewise constant solution. It turns out that 
we must choose this value and place it at * and it must be chosen as the 
state values at a random point between the (See Fig. 2.5). Alter- 
natively (Liu) the values at a point that sweeps out the interval regularly 
will do. 

We assume this set up and state the theorem, show what actually 
happens in some special cases but first a few summary remarks. 

The approach is based on Lax E3l . We are looking for a weak 
solution of 



74 satisfying an entropy condition (stated below). Here u is an n -vector 
and F is a vector valued function. We assume this system is hyperbolic, 
i.e., the matrix F'(u) has n real and distinct eigenvalues for all u in some 
relevant domain. We arrange these eigenvalues, Ayfyi), in increasing or- 
der 



We also assume the system (12.301 1 is genuinely nonlinear in a sense to 
be chosen. A weak solution of (I2.36t means a bounded measurable 
function u such that 



u, + (F(u)) x = 



(2.30) 



A 1 <A 2 <...<A„. 



(2.31) 




(2.32a) 
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for x e an d a weak solution with initial value u(x, 0) = (p(x) means 
a weak solution satisfying 

J {Xtu + X xF(u))dx + J x(x,0)(p(.x)dx = (2.32b) 
f>0 

for all smooth vectors ^ vanishing for large \x\ + 1. A piecewise continu- 
ous solution is a weak piecewise continuous solution and hence satisfies 
the jump condition across a discontinuity: 

S[u K ] - [F k ], k - 1,2,..., n, (2.33) 

where S is the speed of propagation of discontinuity and [•] denotes the 
jump across the discontinuity. 

We now formulate an entropy condition by requiring the following 
to hold: For some k,l < k < n, 

A k {u ( ) > S > A k (u r ) 
while \ (2.34) 

Ak-\(u e ) < S < A k+ i(u r ) 

Here U( and u r are the states to the left and right of the discontinuity 75 
respectively. The eigenvalues A'.s are also called characteristic speeds. 
The condition (I2.34t says the k th characteristic meets the discontinuity 
from the left and the n - (k - 1) characteristic from the right, the total 
being equal to (n + 1) and thus (n + 1) quantities will determine (2n) 
unknowns and the 'shock speed'. This agrees with gas dynamics and 
combustion. 

A discontinuity across which d2.33t and d2.34t hold is called a k-th 
shock and 5 will be called a shock speed. 

Suppose for some k, 1 < k < n, grad l( A k + and is not orthogonal to 
r k the corresponding eigenvector. If this is so, we say k th fan is genuinely 
nonlinear. We normalise r k so that 

r k -gmd u A k = l. (2.35) 

If on the otherhand r k ■ grad I( A k = we say the k th fan is linearly degen- 
erate. 



66 



2. One Dimensional gas dynamics 



We consider an example from gas dynamics. The equations read 
(See ( |2.9(a)| » 

P; + M/7jt + PC U x - 

pu t + p x + puu x = 
S t + S x u = 0. 

76 Here the matrix F'(u) is given by 

' u pc 1 (T 

p- 1 M 

k 

and A - u is an eigenvalue and the corresponding eigenvector is given 

by 

°) 
o 

where uj + is arbitrary. Thus the characteristic field corresponding 
to this eigenvalue is linearly degenerate. (This actually leads to special 
difficulties in computation, see, e.g., Harten [ 17]). 
We now state the main result. 

Theorem 1. The set of states u r which are connected to U{ for \u r — 
U(\ sufficiently small through a k - shock from a smooth one parameter 
family u r = u{e), —e Q < e < 0, u(0) = uf, the shock speed is also a 
smooth function of e. 

Remark. The entropy condition gives the one sided interval. 

We now turn to an important class of solutions, centred rarefaction 
waves; these are the solutions which depend only on the ratio (x—x )/(t- 
t ), x , t a are the centre of the wave. 

Let u be a rarefaction wave centred at the origin: 

u(x, t) = h(xjt). (2.36) 



77 Substituting this in d2.36t . we see that 
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x 1 
— h' + -F'{u)ti = 
t t 



(2.37) 



where ' denotes differentiation with respect to £ = x/t. Thus £ = A^is 
an eigenvalue of F'(u) and h' is a corresponding eigenvector; h is called 
a & - rarefaction wave. 



in view of (12.35b . we can take 

ti = r(/i). 



(2.38) 



Put A - A{u{); ( 12.381 ) has a unique solution satisfying the initial condi- 
tion 

h(A) = ut; (2-39) 

h is defined for all £ close enough to A. 

Let e > be such that /j is defined for A + e; write u r - h(A + e). 

We now construct the following piecewise smooth function u(x, t) 
for t > (See Figure 2.6). 



u(x, t) 



ut for x < At 

h(x/t) for At < x < (A + e)t 

u r for (A + e)t, x. 



(2.40) 




Fig. 2.6. 



This function u satisfies the differential equation (I2.36t in each of 
three regions, and is continuous across the lines separating the regions. 
We shall say that in u the states U[ and u r are connected by a centred k - 
rarefaction wave. 
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Theorem 2. There exists to every state Uf a one parameter family of 
states u r = u(e), < e < e Q , connected to U{ by a k - rarefaction wave. 

We now turn to another important problem, Riemann problem, in 
which the initial value given are 



Theorem 3. There always exists a solution (a fan wave) to the Riemann 
problem if\u$ — u\ \ is sufficiently small. 

The proof uses only the implicit function theorem and we refer to 
Lax t2~3\l. 

Norms: What can we expect: u x e L , not L 2 , at best since u has 
jump discontinuities. This suggests looking for a solution u of bounded 
variation. 

We now state a theorem due to Glimm and for the proof we refer to 
Glimm fEl . 

Theorem 4 (Glimm). Let H2.30t be hyperbolic, strictly (genuinely) non- 
linear and F be smooth in a neighbourhood ofu = v,a constant vector. 
Then there is a K < oo and a 6 > with the following property: 
79 If the initial values uix, 0) are given so thaQ 



then there exists a weak solution of \2. 361 for all x, t > with initial 

values u(x, 0) such that 



u = u for x < 
u = u\ for x > 0. 



d x = \\u(.,0)-v\\ oo + T.V.u(.,0) <6, 




llu-vlU <£||w(.,0)-v|L 

T.V. u(.,t) < K(T.V.u(.,0)), t > 



(2.41a) 
(2.41b) 




(2.41c) 



— oo 



2 T.V. = Total variation 
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For a restricted class including gas dynamics and if 



||w(.,0)-v|Ul +T.V,u(.,0)) < 6 



then there exists a solution satisfying \2.41a\ and \2.4lh\ . 

We now describe an approximate method developed by Glim to 
solve any initial value problem u(x, 0) = u (x) when the oscillation of 
u (x) is small. The solution u is obtained as the limit of approximate 
solutions w/j, as h — > 0, which are constructed as follows: 

(I) Uh(x, 0) is a piecewise constant approximation to u (x) 

Uh(x, 0) = nij for jh < x < (j + l)h, j = 0, ±1 (2.42) 

where ray is some kind of mean value of u a (x) in the interval 



(jh, (j + l)h). 

(II) For < t < hj A, Uh(x, t) is the exact solution of d2.30t with ini- 



tial values u(x, 0) given by (12.421) : here A is an upper bound for 80 
2\Ak{u)\. This exact solution is constructed by solving the Rie- 
mann IVP's, 



j = 0, ±1, Since the oscillation of u Q is small, m,_i and ra ; - are 

close and so by Theorem |3 this IVP has a solution consisting of 
constant states separated by shocks or rarefaction waves issuing 
from the points x = jh, t - (See Figure 2.7). As long as 



these waves do not intersect each other and so the solutions of the 
IVP d2.43t can be combined into a single exact solution u/,. 




(2.43) 



t < h/ A 



(2.44) 
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Fig. 2.7. 



(Ill) We repeat the process, with t = h/A as new initial time in place of 



It is not at all obvious that this process yields an approximate 



show that the oscillation of u(x, nh) remains small, uniformly for 
n — 1,2, . . . and so that one can solve Riemann IVPs (12.43b . This 
estimate turns out to depend very sensitively on the kind of av- 
erage used to compute the mean values mj. Glimm has used the 
following method to compute mf. 

A sequence of random numbers a\,ai, . . . uniformly distributed 
in [0, 1] is chosen; m" then mean value of u(x, nh/A) over the in- 
terval (jh, (j + \)h) is taken to be 



Glimm proves the following. 

Theorem 5. A subsequence of Uh converges in L 1 with respect to x, to 
a weak solution of \2. 3(M . uniformly in t, and for almost all choices of 
{a n }. 

For the proof we refer to Glimm fi"5ll and we illustrate here by an 
example; see Lax [23]. 

Consider a Riemann IVP 



t = 0. 
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solution Ufi which is defined for all ?; to prove this one must 



m n - = u(jh + a n h, nh/A). 



(2.45) 




U( for x < 
u r for x > 
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where U( and u r are so chosen that the exact solution u consists of the 
two states U(, u r separated by a shock, 

(Up if x < st, 
(2.46) 
u r if x > st, 

where s is the shock speed. We may take A > \s\. Assume s > 0. 82 
Glimm's recipe gives 



Uh(x,h/A) = 

where 



\U( if x < J\h 
\u r if J\h < x 



1 1 if a\ < si A 
1 if si A < a\ 
Repeating this procedure n times, we obtain 



Uh(x, nh/A) 



I Uf for x < J n h, 
I u r for J n h < x, 



where J n = number of a'-s, j = 1,2, ... ,n, such that aj, si A. Since {aj} 
is a uniformly distributed random sequence in [0, 1] 

J n s 
n A 

with probability 1 ; this tells the approximate solution tends almost surely 
to the exact solution. 

One would like to prove more about the solution of the initial value 
problem. In particular the uniqueness of the solution. Some results are 
contained in DiPerna [ 8 1 but they are not applicable to the general initial 
value problem because they do not admit shock formations within the 
flow. 



2.10 Combustion. Detonations and deflagrations 

In this section, we present a brief account of the elementary theory of 83 
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detonation and deflagration waves. These differ from shocks as the in- 
creased pressure releases energy and converts one gas into another in a 
chemical reaction. We denote this energy, per unit mass, the energy of 
formation, by g and the total energy 

E = e + g 

where e is internal energy. In this case the equations of the gas dynamic 
are 

p t + (pu) x = 
(pu) t + (pu 2 + p) x = 0, (2.47) 
E t + ((E + p)u) x = 0, 

where p is the density, u is the velocity, p is the gas pressure. But E 
depends on the precise nature of the gas and is given by 

E =pE + pu 2 /2. 

Let the subscript refer to unburnt gas and 1 to burnt gas. Let the 
unburnt gas be to the right of the reaction zone and let U be the velocity 
of the reaction zone. Then the two laws of conservation of mass and 
momentum are identical with the corresponding laws for shock fronts 
and we have 

PoV = Pivi - m, (2.48) 
Po + Povl = pi + pw\ (2.49) 

84 where v, = Uj - U, i - 0, 1, are the relative velocities. The law of 
conservation of energy now takes the form: 

E (o \t , Po ) + p T + v 2 J2 - £ (1) (ri,Pi) + PiTi + vj/2 (2.50) 

where E^ and E^ are two different energy functions. As in the case of 
shock fronts we consider the Hugoniot function 



H(tuPi;t , Po ) = £ (1) (ri,pi) - E {o \t , Po ) + (n - t )( Pi + Po )/2. 
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It should be stressed that the conservation of energy ( 12.501 ) is entirely 
different from that in the case of shock fronts but that any law derived 
from the conservation of mass and momentum still holds. 
Using (|2~3B1 and d2~391 ). ( H30l can be written as 

H(tuPi;t o , Po ) = 0. (2.51) 

As in the previous case, the graph of (ji,pi) in the (t, p) plane, which 
satisfies (12.511) for fixed (r , p ) will be called the Hugoniot curve with 
center (t c , p Q ). For polytropic gases, we have 

e = -, y > 1. 

y~ 1 

If we set A = g - (A < for an exothermic process) and /i = , 

y + 1 

we find 

= 2fi 2 H = - Po (t o -p 2 t0 + pi(Ti -/z 2 t ) - 2// 2 A. (2.52) 

As in the previous case, this is a rectangular hyperbola. This time the 85 
point (r , p ) does not lie on the hyperbola because of the extra term due 
to A. In fact, if A < 0, (T ,p ) lies below the Hugoniot curve; see figure 
2.8. We assume this in general. 

The lines through (t , p ) tangent to H = are called Rayleigh lines. 
Their points of tangency, S \ and S 2, are called Chapman- J ouguet ( CH) 
points. The portion p > p a , t > r , of the Hugoniot curve is omitted 
because it corresponds to the impossible case in which m 2 < 0. Here we 
have used the relation 

2 Po-Pl 
—m = , 

To ~ Tl 

which follows from <!2.48b and A2.49I ). The upper portion of the Hugo- 
niot curve corresponds to detonations (increase in pressure); the portion 
above S 1 corresponds to strong detonations and below S \ to weak det- 
onations. The lower portion of the Hugoniot curve corresponds to de- 
flagrations (decreases in pressure). These two pieces might appear to be 
playing the role of a "state behind" curve and we might anticipate that 
shocks and detonations are linked while deflagrations and simple waves 
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are related. However, considerations of the internal mechanism appears 
to eliminate deflagrations and weak detonations except in special cases. 

The relative speeds of fronts governed by the Hugoniot curve can be 
determined by differentiation. 




CJ point 



Fig. 2.8. The Hugoniot curve for exothermic gas flow 

86 

First of all, 

- dH(r,p) = TdS + -{(r - r )dp + {p- Po )dr}. 

Thus for a Chapman-Jouguet process where (t - T„)dp + (p- p )dT - 
we have dS - 0. So at this point 

P- Po _ dy_ _ dp | i= _J_ c 2 
t - t dr dr Tj 1 
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Thus by the mechanical conditions 

V l = - T l(P ~ Po)/( T - T o) = c\. 

That is, the state behind is sonic relative to the front. On the other hand 

v l = ~ T l(P ~ Po)/(r - t ) 
and along the Hugoniot curve 

dv \ = T l((p - Po)dT - (r - T Q )dp) 

so that v has stationary value for a Chapman-Jouguet process and fur- 87 
thermore it is a minimum. Hence the flow ahead of a strong or weak 
detonation is supersonic and of a strong or weak deflagration is sub- 
sonic provided the Hugoniot curve has the usual convexity. 

By looking at the Hugoniot curve for the state behind in a similar 
way we find that v\ is monotonic on each branch. Hence the flow is 
supersonic behind a weak detonation, subsonic hebind a strong deto- 
nation, supersonic behind a strong deflagration and subsonic behind a 
weak deflagration. See Fig. 2.9. 
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Weak deflogration Strong deflogration 

(b) 

Fig. 2.9. 



88 This leads to Figure 2.9 which shows the way the three character- 

istics can leave and enter a front. Each characteristic is carrying one 
datum. Thus for a strong detonation, one piece of data must be given 
along with the state in front and the two remaining quantities are then 
determined. See the corresponding shock wave problem. 



2.11 Riemann problem with detonations and defla- 
grations 

We assume there are only two gases, burnt and unburnt, and that the 
unburnt gas is on the right with ur - 0. Then there are the following 
possibilities: 
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To the right of the contact discontinuity or slip line either 

a) a strong detonation or 

b) a CJ detonation followed by a rarefaction wave or 

c) a rarefaction wave 

and on the left either a shock or a rarefaction. 

This includes the unlikely case (c) when the pressure in the unburnt 
gas pr exceeds that in the burnt pi. Weak detonations and all deflagra- 
tions have been eliminated for other reasons to be discussed later. 

The S B curve for the left state is the same in section 12.81 The S B 
curve for the right state row consists of the strong detonation branch 
connected at the Chapman- Jouguet point to an adiabatic curve p - p(t) 
and there is also the rarefaction curve through pr for use in case (c). 

Exactly the same kind of argument then shows the Riemann prob- 89 
lem for the case pr < pl can be used again, i.e., connecting the two 
continuous curves by a line segment and matching up the values of u 
which are given by the same formulas unless there is a C - / detonation 
when 

u R = c(t*,S*) - £ C j(t,S*) + £ C j(t*,S*), 
and ul as before. 

Here, t* , S* are the specific volume and entropy at the Chapman- 
Jouguet detonation and £qj is the T corresponding to that state. 

The originale shock wave argument continues to work for pr > pi 
if a shock or rarefaction on the left and a rarefaction on the right are 
generated. The excluded case is the two-shock case because the unburnt 
gas is detonated. For 

u L < [(PL ~ p(f(T R )))(f(r R ) - t l )] 1/2 

there will be a solution without detonations. If ul does not satisfy this 
inequality, then there is a detonation. However, uniqueness is lacking. 
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2. One Dimensional gas dynamics 



2.12 Internal mechanism 

We now investigate why only certain processes take place by looking 
into the internal mechanism of the front. We look for steady state solu- 
tions in one dimension where we now assume the flow has viscosity and 
heat conductivity. Let the temperature be 9. The conservation of mass 
remains the same as before: 

pv = p v = m, a constant. 

90 We seek a continuous flow that moves from a constant state at x = -oo 
to a constant state at x = +oo. We assume the flow is moving from 
right (unburnt) to left (burnt) and the reaction front has fixed velocity, 
so m > 0, and x = +oo is ahead, x = -oo is behind. The conservation of 
momentum is given by 

2 dv 
pv + p- u— - p, a constant, 
dx 

where p. is the coefficient of viscosity. 

Next, at any stage the gas has internal energy E^ £ \6) and the energy 
balance is influenced by heat conduction. Denoting by A, the coefficient 
of heat conduction, the energy balance is written as: 

dO r c \ 1 7 dv 

-A— + m{E '(9) + -v } + v{p - p—\ - mQ = constant. 

dx 2 dx 

Following Friedrichs (see [11]), we assume the balance between 
burnt and unburnt gases is 

-v^+(l-e)S(0) = O, 
dx 

where we try to avoid specifying too much about S(9). So, we have 
three autonomous equations and they have singular points exactly at 
possible end states. 

We write E (0 \9 ) = E and E (l \9\) = E\ and investigate the sin- 
gular points and find that at x - -oo (state '0') there are in general 
solutions behaving like 



e aix^ e a 2 x^ e a 3 x^ 



2.12. Internal mechanism 
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where ce\ > 0, 02 - and 

a 3 > if v > c„ 
a 3 < if v c < c . 

Here the suffix corresponds to the unburnt gas and c is the sound 
speed in the unburnt gas. 

The manifold of regular solutions (i.e., tending to constant as x — > 
-00) has two free parameters if v > c and only one if v < c a . Call 
the number of parameters y . Analogously, at the other end, we find the 
manifold of regular solutions has one free parameter if v\ > c\ and two 
if vi < c\ ; call the number of free parameters y\ . 

When can we even hope to find a solution going from state (0) to 
state (1)? We have three quantities to solve for. We must impose (3 —y ) 
initial conditions to leave the state (0) and y\ parameters characterise 
regular solution at +00. But one is used up by an arbitrary shift in x. 
Hence (71 - 1) parameters are to be chosen subject to (3-y ) conditions. 
So we need y\ - 1 > 3 - y a , i.e., 4 - y - y\ < or else some other 
quantity must be specially chosen. Thus: 

Strong detonation v a > c„, y () - 2 

4 - y„ - 71 = 

vi < ci, 71 = 2 
Weak detonation v > c , y„ = 2 

4 - y„ - 71 = 1 

vi > a, 71 = l 
Strong deflagration v < c , y - 1 

4 - y ~ 71 =2 

vi > ci, 71 = 1 
Weak deflagration v < c , 7 = 1 

4 - 70 - n = 1 

vi < ci, 71 = 2 
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2. One Dimensional gas dynamics 



So over determinacy the internal mechanism matches the under de- 
terminacy of the characteristic problem as given in § 12.1 II 

When the topology of the solution curves is worked out, it turns out 

that 

(i) All strong detonations are possible, 

(ii) CJ detonations are only possible if one of the parameters satisfies 
a special condition, 

(iii) Weak detonations are not possible except under very special con- 
ditions, 

(iv) Weak deflagrations are as in (iii) and 

(v) Strong deflagrations do not exist. 



Chapter 3 

Two dimensional steady flow 



3.1 Equations of motion 

The next in simplicity to one dimensional flow is a two dimensional 93 
steady flow which is also irrotational. We use the following notations 
throughout this chapter. 

Let p be the density of the gas and p be the pressure. They are 
considered as functions of the cartesian co-ordinates x,y. Let u, v denote 
the velocity components along the ^-axis and y-axis respectively. The 
equations of motion reduce to 

Conservation of mass: 

(pu) x + (pv) y = 0. (3.1) 

Conservation of momentum: 

(pu 2 ) x + (puv)y + p x - (3.2a) 
(puv) x + (pv 2 ) y + p y = 0. (3.2b) 

We restrict ourselves to situations with weak shocks and assume 
the flow is isentropic, i.e., p - p{p) with p'(p) > 0. The irrotational 
condition implies 

u y = v x . (3.3) 
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3. Two dimensional steady How 



With these assumptions the equations for momentum reduce to 

1 2 

UU X + VUy H C p x = 0, 

P 

1 2 

UV X + Wy H C Py = 0, 

P 

94 where c = {dp I dp) 1 ^ 2 is the sound speed. Equivalently, 

V (I( M 2 +v 2 ) + J jdp) = 0, 
which in turn implies Bernoulli's law 

X -cf = X -(u 2 + v 2 ) + i(p) 
c 2 

is constant, where i(p) = f — dp and q* is the speed at zero density. 

p 

We recall that across a steady shock the following relations hold: 

[p{3 - l)) ■ ft] = 0, (3.4) 
[(3 - U) x ft] - 0, (3.5) 
- U) ■ ft) 2 + p] = 0, (3.6) 

d\4-U\ 2 + e(p,p) + ?-] -0, (3.7) 
2 p 

where - («, v) is the velocity vector, u is the shock speed, ft is normal 
to the shock and [•] denotes the jump across the shock. Using the fact 
that the entropy and vorticity changes are to be neglected, see Chapter 
EJ we see that (I3.7t may be written as 

[\\q-U\ 2 + i(p)] = 

or equivalently, 

[q* 2 ] - [fl • H = 0. 



3.1. Equations of motion 
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In a frame of reference where the normal component of the shock ve- 
locity is zero, we have the tangential component of q is continuous and 
so [q] ■ U = 0. Hence q* is continuous. This, along with the assump- 
tion that the entropy is constant, replaces the conservation of energy and 
the normal momentum equation. In fact energy and normal momentum 
equations are not conserved. We are left with the two shock conditions 
given by the conservation of mass and the continuity of the tangential 
component of velocity. 

Thus we need to consider the flow which satisfies the conservation 
laws: 

(pu) x + (pv) y = 0, (3.1a) 
Uy - Vx = 0, (3.3a) 

in their integrated form where p as a functions of the speed q is given by 
Bernoulli's law 

1 1 2 

-q 2 + i(p) - 2^*"' wnere q 2 - u 2 + v 2 . 

Equation d3.3al ) implies there is a functions (f>(x,y), called the potential 
function, such that 

<f> x = U, (fry = V. 

Similarly, (13.1 at implies, there exists a function \fr(x,y), called the stream 
function, such that 

iffy - PU, ~l/f x - PV. 

The shock conditions reduce to the conditions: 

<p, \\i are continuous. 

From the definitions of <f>, if/ we see that if we introduce w = u - 
the equations of motion reduce to 

V^dil/ 

d(f> + „ - wdz, (3.8) 

p(M) 

where z = x + y—ly, and p(\w\) is given by Bernoulli's law. From d3- 8I > 
any number of alternative equations are easily written down on the basis 
of the perfect differential properties of (13.81) as we shall see in the section 
on hodograph transformations. 
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3. Two dimensional steady How 



3.2 Classifications of flow equations 

From Bernoulli's law, 

c 2 

(u du + v dv) H dp — 0. 

P 

Using the mass equation (13.11) and the above we obtain the equation 

(c 2 - u 2 )u x - uv(u y + v x ) + (c 2 - v 2 )v x - (3.9) 
or for the potential <p the well-known governing nonlinear equation: 

(C 2 - U 2 )4> xx - lUV^xy + (C 2 - V 2 )<Pyy = 0, (3. 10) 

where c, u, v are functions of V0. It is convenient to introduce the Mach 
number M: 

M = q/c. 

The characteristics of the partial differential equation ( 13.101 are given 

by 

dy _ uv ± c 2 VM 2 - 1 (3 11) 

dx c 2 - u 2 

97 The equation A3. 101) is hyperbolic, elliptic or parabolic for M > 1, M < 1 
or M = 1 respectively. In the first case, the flow is said to be supersonic, 
in the second subsonic and in the third sonic. Among these, the first 
two types of flows are more or less thoroughly studied and the theory 
is understood if not complete. But when both M > 1 and M < 1 occur 
in a single flow, we call it transonic. This mixed case has many open 
problems. 



Remark. By Bernoulli's law, 



c 2 



qdq H dp = 0. 

P 



Therefore, 

d q 2 
-j-ipq) = p(l - -j) 

dq c z 

At the sonic line q 2 = c 2 , therefore, as in the case of a scalar conserva- 
tion equation (Chapterffl, pq has a maximum. 



3.3. Supersonic Flow 
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3.3 Supersonic Flow 

The supersonic case is equivalent to the hyperbolic systems we have 
already studied where we may treat x or y as a "time" variable. Thus 
Cauchy problems may be solved. The characteristics are given by d3. 1 lb : 

dy _ uv ± c 2 VM 2 - 1 
dx c 2 - u 2 

By looking at the case v = 0, the path of a particle is given by 
dy/dx = 0, we see that the path bisects the characteristics. Note that as 98 
M — > 1 the characteristics become perpendicular to the direction of the 
flow and tangent to each other. 

There are simple waves, Riemann invariants and a solution to the 
analogue of the Riemann problem and the piston problem. However, 
the two elementary flows of greatest interest correspond to the flow past 
a bend in the wall (see the figures below). 




Fig. 3.1. 
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3. Two dimensional steady flow 



A continuous flow is described by a simple wave whenever it is ad- 
jacent to constant state. An inward bend (Fig. 3.1(a)) causes the char- 
acteristics to form a cusped envelope and hence a shock. An outward 
bend (Fig. 3.1(b)) has a continuous rarefaction wave possibly ending at 
99 zero density and the escape velocity which, by Bernoulli's law, is equal 
to q* . A sharp straight bend yields a rarefaction wave or a straight shock 
(Figures 3.2 (a) and 3.2(b)). 



These problems are easily solved by noting that the image of the 
rarefaction wave if a characteristic in the (u, v)-plane and the image of 
the possible states behind the shock is given by a so called shock polar. 
It is convenient to look at these in (9, q) - plane, where 9 is the flow angle 
tan -1 (v/w). 

3.4 Shock polar 

From the discontinuity conditions in the form 

[pu]dy/dx - [pv] = 0, 
[u] + [v]dy/dx = 0, 

100 where dy/dx is the slope of the shock, we find that with an initial state 
(q a , 0) with p - p(q ), if the state behind the shock is (q cos 9, q sin 9), 
then 




(a) 



(b) 



Fig. 3.2. 



q cos 9 = 



pg 2 + p () ql 
q (P+Po) 



(3.12) 



3.4. Shock polar 
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In the limit of a weak shock, we note that a shock can become weak 
in two ways: either the shock becomes characteristic or the flow on both 
sides becomes sonic. To see this, we set 

F = pq, F = p„q , F/F Q = 1 + 6G and q/q = 1 + Sp. 

Then we obtain from A3. 121) 

5p ■ 5G 



cos 9 = 1 + 



2 + dp + 8G 



or 



, e i SpSG 

sin - — 



2 2 2 + Sp + SG 
To a first order approximation, we have 

SG = Sp ■ p • dF/dq. 

Hence, if the shock is weak we have # — > and dp — > 0, and we obtain 



g 2 = (6p) 2 (-dF/dq), 



and the shock is characteristic with 

9= ±6p(-dF/dq) l/2 , 

provided dF/dq + 0. However, if the velocities on both sides of the 
shock are close to sonic, if we approach the limit appropriately: 

9 2 = {5pf. 



Note that in the normal shock case 9 = and 6G = 0. The shock polar 
for the polytropic case is as illustrated in Fig. 3.3 (See Bers [2]). 
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3. Two dimensional steady flow 



n 




1" 



Fig. 3.3. SHOCK POLAR 



An important problem is the detached problem. Suppose a projectile 
or wing is moving in a fluid with supersonic speed. If the projectile (or 
wing) is round the speed vanishes at the tip (stagnation point). There is a 
shock in front, not an attached shock. The shock is curved and the flow 
is constant in front because of the curvature behind the shock, the flow 
is not strictly isentropic or irrotational. If the shock strength is moderate 
a situation that occurs if the speed at infinity is not too high we may 
still use the conditions of irrotationality and isentropy. Then there is a 
subsonic region around the nose. 



Mi > 1 
► 



Shock 




Fig. 3.4. DETACHED SHOCK OVER A BLUFF BODY 



3.5. Equations in the hodograph plane 
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The problem is to solve d3.lt and ( 13.3b . a mixed system, with bound- 
ary condition udy - vdx - on the object. The shock is a free boundary 
in the flow represented in the u, v place by the shock polar. From <I3. 12l > 
one sees that the flow behind the shock is subsonic on the axis, so there 
is a region of subsonic flow. However, like the nozzle flow discussed in 
§12 there are no transonic difficulties such as those described in §7-11. 
See Bauer et al. Q. 



3.5 Equations in the hodograph plane 

By the hodograph plane, we shall mean either (u, v) or (q, 0) or even 
(cr(q), 6) plane, whichever is convenient since they correspond to simple 
mappings except near q - 0. However near q = 0, the equation for the 103 
potential cp behaves like = and we have a very complete under- 
standing of the flow. The special function cr(q) is chosen to simplify the 
equations. From (13.81 1. we obtain 

P 

Hence the left hand side is a perfect differential. Considering <f>, if/ as 
functions of <x, 6, where cr = cr(q), we find 

q- l e^\<p e + ^^)dG + q~'e^\K + ^^)dcr 
P P 

is a perfect differential. This requires 



Thus, 



(g-'x^+p-y y^o, (3.i3) 

(p"V)«r^ - ?~V«r = 0. (3.14) 
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3. Two dimensional steady How 



For the transonic range, it is convenient to introduce 

9 

cr - J" q~ l p dq. 
Then equations (13. 13l > and (13. 14l > reduce to 

fa = ft<T, 

and 

(Per = -K(o-)lff d , 

where 

1 d(pq) 
p i q 1 dq 

is a function of cr only. Thus 

K(o-) fgg + i/frcr = 0. (3.15) 

Note that the equation is elliptic or hyperbolic according as K(o~) is 
positive or negative, i.e., as cr is positive or negative. Note also that the 
characteristics of (13.151 are given by 

cr 

9 = ± j' yJ-K(cr)do- + constant. (3.16) 



The Legendre Transformation: In solving perturbation problems, 
it will be convenient to use the Legendre transformation which we in- 
troduce now. Recall equations d3.1at and d3.9t . If we regard x,y as 
functions of u, v the equations reduce to 

x v -y u = (3.17) 

(c 2 - u 2 )y v + uv(x v + y u ) + (c 2 - v 2 )x u = 0. 

By ( 13.171 . there is a function x(u, v), Legendre transformation, satisfy- 
ing 

x=Xu, y =Xv 



3.6. Small Disturbalce equation 
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The relation between the potential ft and x is then given by 

cf> = xu+yv-x 

as is easily seen. Similarly, we can introduce a Legendre transformation 
X such that the stream function ft can be written as 

ft = -pxv+pyu -x- 

It is clear that we can consider both^-,^ as functions of q, 9. 
Combining the definitions ofx,X with P.8l l. we find that 

dx = xdu + ydv = x(cos Gdq - q sin 9d0) + y(- sin 6dq - q cos 9d6), 

d(pq) 

dx - yd(pu) - xd(pv) = y(cos 9 dq - sin 6 ■ pqdff) 

dq 

d(pq) 

- x(- sin 9 dq - cos 9pq d9). 

dq 

From which it follows that 

d(pq) i i 
Xe = - (-^-) Xq = (P<V Xe- 

We set 

d& = dqlpq and K = pq—^r~ 

dq 

and the equations become 

Kxe = -%&, X& =Xo, 

where K, a are not the same as K and o~ but have the same basic prop- 
erties. 

Furthermore, 

&Xee+XM = 0- (3-18) 

We shall use this equation later in connection with perturbation prob- 
lems. 

See Bers [2] for more details on these equations and for early work 
on subsonic and transonic theory. 
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3.6 Small Disturbalce equation 

Von Karman's model nonlinear mixed equation 106 

<Px<f>xx + (f>yy = 0, (3.19) 

elliptic for <p x > 0, hyperbolic for <p x < 0, can be derived from ( 13.81 ) by 
expanding 9 about '0' and q about the sonic value c () . From 

pdtp + V-L# = pqe~ ^~~ ie dz 

we have with <f> = c x -<p,if/ = p c ,c y + if/, 

(Po+p- Po) (c dx - d(p) + i(p c dy + difr) 
1 d 2 (pq) 2 

- (p q + t - 9 \o (q - c ) + ■ ■ .) (l - 10 + . . .)dz 

l dq z 

so that from the highest order terms 

-p d(f> + idifs = [A(q - c ) 2 - p o q o i0]dz - B(q - c Q )dx, 

where 

A = -d 2 (pq)/dq 2 \ and B - dp/dq\ ■ c B . 
Here, we have used 

d(pq)/dq\ = and q G = c D . 

Thus 

p <p x = B(q - c ) - A(q - c ) 2 

<Px = -p q Q 
-Po<P y = p q Q 

iffy = A(q - c ) 2 . 

107 The small disturbance equations are obtained by eliminating 9, q and 
thus 

<Px = Po 2 B 2 A~ l fty, to first order in ifr y cc (q - c ) 2 , 



3.7. Transonic How 
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Po<Py = fix- 
But A < so that after rescaling (x -> je and j -> (-|p~ 1 .B 2 A~ 1 )~ 1/ ' 2 ;y) 
we obtain 

02 = { -L p -^B 2 A- l ) l/2 (-2p- l )fi y , 

Po (-^B 2 A~ l )- l/2 <p y = fi x . 

These equations in turn yield d3.19t . 

On the other hand, by expanding the hodograph equations we obtain 

crfia - -fie 
fie = fio-- 

3.7 Transonic flow 

We limit our discussion to a few problems centered around transonic 
wing flow. 

The first question is whether we can find a wing shape which at a 
prescribed subsonic speed at infinity has a smooth flow. Experimental 
observation suggested in the forties that perhaps such a steady flow did 
not exist as the speed came close to sonic at infinity and locally super- 
sonic in a region next to wing. However, Lighthill 031 showed how to 
construct such a wing shape. But this wing was not constructed possibly 108 
because the prevailing sense was that in any case it would be unstable. 
Frankl' and Guderley, [10], [ 16 1, proposed that the explanation for the 
instability lay in the fact that the boundary value problem was ill-posed 
in the sense of Hadamard if the flow was required to be smooth. This 
we proved by the author [29]. 

The implication was that in general flows with transonic regions 
would have shocks, not detached as in supersonic flow but arising in 
the supersonic region or cutting it off. This is in fact the case and the 
shocks have a strong effect on the drag. However, this effect is much less 
than the drag produced at supersonic speeds by the datached shocks. So 
that it has proved very useful to use such wings. In 1331 . Nieuwland has 
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designed an algorithm for finding, not a transonic wing, but a smooth 
transonic cross-sections. 

In the following sections, we shall discuss the relevant boundary 
value problems and the related general theory of mixed equations. Then 
we shall describe the design of aerofoil shapes developed by Garabedian 
et al., and the method introduced by Murman and Cole for finding the 
flows at off design Mach numbers. 

3.8 General theory of boundary value problems for 
mixed equations 

The mixed equations were first investigated by Tricomi, see [36], for the 
equation that bears his name and is the hodograph equation for the small 
disturbance equation: 

Olffgg + tfw = 0, (3.20) 

see the end of § 7. 

A sample theorem on a boundary value problem which illustrates 
that the standard Dirichlet problem would be over determined is: 



o 

u 
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Fig. 3.5. TRICOMI PROBLEM 



Theorem. Suppose i]j satisfies A3. 20t . where i{/ is prescribed on the curve 
C3 and the characteristic C<i. Let D be the region enclosed by the curves 



3.8. General theory of boundary value.. 
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C3, C[ and C2 (C\ is also a characteristic). Suppose the are C3 is star- 
lik^ with respect to the origin: 



Under these conditions, the solution is unique for if/a-, \jjg continuous 
throughout the closure ofD (See Garabedian M 3V ). 

It is reasonable to expect such a theorem. Consider the still simpler 
mixed equation (Lavrente'v, Bitsadze l2"2"l). 



and the same boundary data. Here C\ and C2 are the straight character- 
istics. Let F be the value of if/ on cr = 0. Then G = dip/dcr\cr = found 110 
by solving an elliptic problem is a linear functional of F, in cr ^ 0. 
On the other hand solving the wave equation, for cr < 0, one finds that 
the data on C2 alone determine another linear relation between F and 
G. It is reasonable to expect that one can eliminate G and solve for F 
uniquely. Then tp is easily seen to be unique. 

There are a variety of methods for proving the theorem and also 
establishing the existence of the solution. The method we shall use is 
by an estimate. Without loss of generality, we can consider the nonho- 
mogeneous equation with homogeneous boundary data. We rewrite the 
equation as a system by setting 



9dcr - o-d9 > 0. 



(sgn o-)ifj ee + tfro-a- = 



(3.21) 



on = (a>i,a>2) with a>\ = i//q, 0J2 = ft, 



and thus 



Acog + Bojfj- = f 

(jO\d6 + (jj-idcr = on C2 + C3 



(3.22) 



where 



A = 



0- 
-1 



) 



B - 



1 

1 



J 



(3.23) 



'This condition means, as a point moves along C3 inthe counterclockwise direction 
the angle which it makes the 0-axis is increasing. 
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which is a symmetric system. 

We seek a matrix C with a certain property: Take the 'scalar product' 
of (13.22b with C, where AC and BC are symmetric. Integrate by parts 
over the domain D. Apply the boundary condition. Choose C so that 
the boundary and area integrals are positive definite. Denote the area 
integral by Q(co, oS). Then 

Q(co, to) {f, Ceo), 
where (, ) is defined for any two column vector functions 

f^iflJlY, g = (gl,g2)' 

by 

(f,g) = ff (/igl + hgi)dx dy. 

D 

Thus if / = then co = which proves uniqueness. 
We now proceed to find C. If 



C21 C 2 2 



we find 



1 1 

(Caj,A(jJs) - (ACco,coe) = -(ACo>, u) - -((AC)eOJ,(o), 

provided AC is symmetric, i.e., crcn = -c%\, and 

1 1 

(Ceo, BtOo-) - -(BCoj, - -((BC))crco, oj), 

provided BC is symmetric, i.e. c\\ = C22- Hence the boundary term is, 
by Green's theorem, 

1 r 2 2 

- I co^ibder + cdG) +2u)\U)2(crcdo- - bd6) - 0J 2 (bdcr - cd6). 



Here c\\ - C22 - b and c\2 - c. 



3.8. General theory of boundary value.. 
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The contribution to this term from C 2 + C3 where \p = and we may 
write 

co\ = adcr, a> 2 = -adO (3.25) 



112 is 



2 (bdcr - cdff) {ado 2 + d9 2 ). 



C2+C3 

And the integral on the characteristic C\ where d6 + -^-crdcr = is 



crtoi - OJ2) 2 (b + c ^-cr)}do-. 



Ci 

Hence to fulfil the required positivity C must satisfy: 
(AC)g + (BC)o- is a negative definite matrix 

bdcr -cd6>0 on C 2 + C 3 
(b + c V-o")^o" < on C\ . 

We then find the explicit requirements, 



crbo - (crc)o- ^ 

-b e + Co- ^ (3.26) 



(o-c e + baf ^ {(rbg - (crc^X-bg + Co-) 

in D and the boundary conditions 

(bdcr - cdff) (0-do- 2 + dG 2 ) ^ on C 2 + C 3 
(b - V^c) ^ on Ci 

The choice ft = 9, c = cr for cr > 0, c = for cr ^ satisfies these 
requirements. This completes the uniqueness theorem. Note that the 
requirements would also be fulfilled if C 2 was not characteristic but sat- 
isfied only crdcr 2 + d9 2 > and that the star-shaped condition on C 2 113 
could be freed up by changing b and c. 

A weak solution a> € H t of (13.22b and d3.25t satisfies 

-{Lv,co) = {v,f) 
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for all smooth vectors v such that 



vi = on C 2 + C3 



and, since the matrix A + ^f^aB is singular, such that 



V2 - v\ V-o" - on C] , 



where L = A 3/ 39 + Bd/dcr and //* is defined below. 

This is the adjoint problem. To use the projection theorem it would 
suffice to find a Hilbert space for which (v, f) would be a bounded linear 
functional of Lv for v satisfying boundary conditions. But (v,f) is a 
bounded linear functional in v in some weighted L 2 -space and / in its 
corresponding space. So, we have to find an appropriate space in which 
v is bounded in terms of Lv if 



We proceed to the details. 

Let H* be the Hilbert space of all pairs of measurable functions u - 
{u\, ui) for which the norm 



114 where r 2 = 6 2 + 2 . 

Let V be the set of all function cd - (a»i, 0J2) with continuous deriva- 
tives and such that 



vi = on C\ + C 2 + C3. 



(3.27) 




D 



is finite; the inner product is given by 




D 



co - (0, 0) at r - 0, 
o>i = on Ci + C 2 + C 3 , 



3.8. General theory of boundary value.. 
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and 

IT (-(Lo))l + (Lto) 2 2 )d9 do- < oo 

D 

Let H* denote the Hilbert space of all measurable functions u — (u\, u<i) 
for which 

Hull* = { ff {- u \ + u\)de do-} 112 

D 

is finite; the inner product is 

(w, v)* - J"J*(~ UlVl + u 2 v i)d0 dcr. 

D 

Note that LV c H* . We now state the following. 

Theorem . There exists a weak solution of <\3.22t and ( 13.251 1 for every 
feH*. 

To prove the theorem we require the following lemma which will be 
proved later. 

Lemma. For all v e V, f e H* 

\(v,f)\SB\\Lv\\*\\f\\*, 

where B is a constant. 

Proof of the Theorem; For v e V define 

G(Lv) = (v,f). 

By the lemma G is bounded on LV c H* . Thus by Hahn-Banach the- 
orem G can be extended to H* as a bounded linear functional. Thus by 
classical Riesz's representation theorem there is a t e H* such that 

(Lv,t) = (v,f) for all v e V. 

The function a> defined by <x>i = -t\l£l, a»2 = -t% will belong to H t and 
satisfies -(Lv, to) = (v, f). Thus to is the required weak solution and this 
completes the proof. 
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3. Two dimensional steady How 



Proof of the Lemma: By Schwarz's inequality, we obtain 

i(v,/)uiiviy/ir. 

Thus the proof will be complete if we prove the a-priori estimate 

IMI* ^ B\\Lv\\*, 

where B is a constant. 

We proceed to do this, Set 

v - Cv. (3.28) 

Again consider 

(Lv, v) = (Avg + B V<T , v) - (A(Cv)g + B{Cv) a , v). 

Again by rearranging terms properly, we can integrate by parts. The 
boundary condition d3.7t becomes bv\ + CV2 = 0. Set 

(Lv,v) = h+h, 

where l\ is area integral and I2 is surface integral. It is easy to see that 
if (13.261 is satisfied then I\ is positive definite; also I2 is non-negative. 
Thus 

|(Lv,v)|^/i. 
On the other hand, for any A > 0, 

|(Lv,v)|^i||Lv|r 2 + i||v||, 

and hence 

h S A\\Lv\f + j\\v\t A>0. 

Thus if llvll* can be estimated in terms of I\ then by choosing A suf- 
ficiently large, we can estimate ||v||* in terms of ||Lv||*. For the same 
choice of b and c one can estimate ||v||h, in terms of I\. This estimate is 
obtained by the same method as was used in the uniqueness theorem. 
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Once the existence of the weak solution has been established, one 
proceeds to determine whether it has in fact some better properties but 
we shall not do that here except to say that since the smoothness prop- 
erties are local, elliptic methods suffice in elliptic regions, hyperbolic in 
hyperbolic regions, and that the solution is a strong solution everywhere. 

For general description, see Morawetz ODTl . See also Osher l35l for 
a different approach. 

3.9 The boundary value problems of transonic wing 
flow 

The two most important boundary value problems for transonic wing 
flow are for the nonlinear solution and for the linearized flow about it. 

In the physical plane, the flow potential tp satisfies O.lOt . the bound- 117 
ary condition d<p/dn = on the wing and at infinity, V0 is prescribed. 
We know from incompressible flow that we need, for a well-posed prob- 
lem, the additional condition (the Kutta-Joukowski condition) that the 
circulation at infinity adjusts itself so that the flow past a wing with a 
cusp at the trailing edge has finite velocity. For compressible flows by 
Bernoulli's law, no infinite velocity is possible any way and the circula- 
tion adjustment is chosen to prevent. 

It is not unreasonable to anticipate that the problem is overdeter- 
mined on the basis of the boundary value problems discussed earlier. 
There are two possibilities: Shocks in general and special smooth solu- 
tions. 

The problem with shocks for a general aerofoil has only been tackled 
numerically. There is some indication that its perturbation problem is 
well-posed, see Morawetz 1 1271 . ] 

There exist numerical codes ("analysis" codes) for solving the 
boundary value problem using artificial viscosity or penalty methods. 
The original viscosity method of Murman and Cole [ 32 ] is described in 
§ 12. It's basic ideas have been incorporated and considerably modified 
by Jameson [18] and can now be used on three dimensional problems. 
Bristeu et al. [4] has used finite elements and a penalty method. 
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3. Two dimensional steady How 



The lag of the theory behind numerical experiment is not suprising 
especially when one realizes how limited the theory is with respect to 
one dimensional flow. 

In the case of no shocks, the problem can be looked at in the hodo- 
graph plane, under the hodograph transformatioru see § 5, the prob- 
lem transforms into a free boundary value problerro Consider the sym- 
metric wing section which has a singly curved image in the hodograph 
plane. The equation, say for iff, ( 13.151 . is linear. There is a prescribed 
singularity at the image point of the point at infinity. The boundary con- 
dition ip = is imposed on the axis 8 = and the unknown image of the 
wing, see Fig. 3.6. But there is a second boundary condition because 
the flow angle 6 is a prescribed function of the arc length on the wing. 
Using (13.81) . this yields, on the free boundary, 

d(p = qi^e^idXiO) + ^P\dY(6)), 

where x = X(6), y = Y{6) describe the wing. This is only one condition 
since tan G = dY/dX. Then using 

d(f> = tyo-dO - Kif/odcr 

we have the extra condition on \p: 

-ipo-dO + KiJ/gdcr = q(cr). 

Solutions for special shapes known as super-critical airfoils cor- 
responding to smooth flow pasta wing can be found. One solves the 
boundary value problem without the last condition (3.31) using some 
smooth boundary and a well posed problem. But we know from § 8 that 
we cannot expect to solve the problem with full Dirichlet conditions. 
Instead use the following procedure: 



2 This approach has been explored by Brezis and Stampacchia 1 3 1 for subsonic flows 
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Fig. 3.6. 



Solve the boundary value problem: 

Kif/ eg + if/o-o- = in D 

<A = on C1+C2 + C3 + C4 + C 5 , 

with prescribed singularity at (0, cr M ). Here C\ + C2 and C3 + C4 are 
smooth arcs satisfying 

crdo- 2 + d8 2 ^ 0, 

and C5 is a slit on cr-axis s.t. crco S cr < 00 on C5. 

The rest of the boundary of D consists of two characteristics T_ 
and T + issuing from the origin until C2, C3 are intersected. It can be 
shown by the methods of the preceding section that this is a well-posed 
problem. The singularity at (0, (Too) is a bit messy (see Gilbarg [14]; 
but for our purposes it suffices to treat it like the incompressible flow 
singularity which would require taking ^(croo) = 1 and 

(Ae - V^TiAtr - A((f> + V^Tcr) 372 = O(|0 + V^Tcr^ 2 ) (3.32) 

where A is related to the prescribed speed of the wing). 
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3. Two dimensional steady How 



We now have part of a flow and part of a supercritical wing image 
(Ci +C2 + C3 + C4) in the hodograph plane. To choose the wing, we con- 
tinue the solution ip across the characteristic gap, T + and r_, by solving 
the appropriate Goursat problem. It is not unreasonable to expect that 
there will be curve joining C2 to C3 on which i{t = especially if the 
hyperbolic region is small. At this stage we have solved the hodograph 
problem. The next stage is to find the physical image using (I3.8I) . A 
short calculation shows that this will fail if K(cr)i//g + ify changes sign (a 
limiting line occurs where K(o~)if/i + ify - 0). But again for sufficiently 
small hyperbolic regions this does not happen and we do in fact find a 
physical flow. 

In the next section, we describe the method used by Garabedian to 
generate smooth flows. Fung et al. [ 12] has used the idea of finding first 
a purely subsonic flow with a sonic line separating two subsonic regions 
adjusting the equation of state and then continuing the flow from the 
sonic line into the smaller region using the right equation of state and 
finding a new profile. 

There are other possibilities for generating supercritical wing sec- 
tions all involving some form of unique continuation. 

Having constructed a smooth flow and a wing section, one way or 
another, one asks what happens to this flow when it is a disturbed say 
by changing its tilt (angle of attack) or by changing its speed at infinity 
(Mach number). The evidence both numerical and experimental is that 
a shock develops at the rear sonic point on the wing profile. It increases 
the drag (and hence the fuel consumption). Theoretically there are no 
results on the nature of this shock flow. 

3.10 Perturbation boundary value problem 

It would be useful to know what flows close to supercritical (smooth 
transonic) flows look like. First one looks at the general perturbation 
problem assuming that the nearby flow is smooth. This leads to a con- 
tradiction. Then one asks what actually happens. 

First, the perturbation equation has to be determined. The direct 
method from equations ( 13.91 ) is tedious. Instead, for the Legendre po- 



3.10. Perturbation boundary value problem 
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tential 

X(u,v) = xu+yv- <p(x,y) 
we see that the perturbation Legendre potential 8x satisfies 

8x = Sx.u + Sy.v - <p x 8x - <p y 5 y - Sep 

in terms of a perturbation potential 8<p and the perturbations 8x, 8y. But 
cp x - u and <p y = v so that Sx = -8<p to first order. Thus, in the hodograph 
variables 0, & of the undisturbed flow since 8x satisfies (13.181) so does 
8(p. Thus 

K(8<p) ee + (8<p) && = 0. (3.33) 
On the perturbed boundary given by 122 

y = Y(x) + 8Y(x) 

we find 

dd> d 

-f(x, Y + SY) + —{8$) = 0, 

on on 

or to first order 

— - (— <f> y )8Y - 0. (3.34) 
on on 

At infinity, if there is a change in Mach number at infinity, there is a 
prescribed singularity of order 3/2 as in the unperturbed case; see d3.32t . 

If we restrict ourselves to solutions with continuous derivatives (no 
shocks) then one finds by using the methods of the first section of this 
chapter, that the problem is ill-possed if we fix the Mach number and 
change the profile, i.e., 8Y + 0. For the symmetric profile see Morawetz 
[28], for the non-symmetric case Cook [5]. Very partial results exist if 
we change Mach number, Morawetz [ 29| III]. 

The simplest proof amounts to showing that the solution is uniquely 
determined up to a one parameter family by the prescribed data ourside 
a "characteristic gap", i.e., by data on C\ + C2 + C3 + C4, see Fig. 3.6, 
since cr = o"(cr). Therefore the function 8Y(x) cannot be arbitrary in 
the gap since 8<p is determined by unique continuation from Y + , V . A 
more elaborate proof shows that 8Y(x) = in the gap if 8Y = on 

C1+C2 + c 3 + c 4 + c 5 . 
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3. Two dimensional steady How 



The next question is how to find a well-posed perturbation problem 123 
that represents a disturbance with shocks. The answer probably lies in 
finding a suitable singular perturbation for the case where on the whole 
boundary the value of 8Y is given arbitrarily. This might be accom- 
plished by admitting singularities into the perturbation velocities at the 
places where shocks are expected, i.e., the points where the sonic line 
hits the profile. Thus the perturbation flow velocities would be a small 
variation on the unperturbed flow velocities but not right at the sonic 
points on the boundary. See iBTl where such a singular Dirichlet prob- 
lem is solved. 

3.11 Design by the method of complex characteris- 
tics 

The method of complex characteristics has been introduced and used 
successfully by Bauer et at in the computation of flows and pro- 
files. It began with the computation of flows with Mach number greater 
than one at infinity where the object was to determine the subsonic flow 
behind an analytical shock, see § 4. Its full strength came in the com- 
putation of supercritical airfoils, i.e., transonic but shock free with some 
Mach number less than one at infinity. We sketch here the principles 
involved in one of the early computations. 

We are given a velocity at infinity for a flow. The object is to find an 
airfoil and a smooth flow past it with this velocity at infinity and with 
somewhat indefinitely specified characteristics: 

124 (i) a large supersonic region, 

(ii) a large decrease in the pressure (a large increases in the velocity) 
to control boundary layer separation on the forward end, and 

(iii) a subsonic cusped trailing edge with the streamline from the upper 
surface meeting the one from the lower surface smoothly. 

In the most recent work the object has been to specify the speed as 
a function of arc length along the airfoil but we will not describe this. 



3.11. Design by the method of complex characteristics 
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We write the equations of motion ( 13.31 ) and ( I3.9t as 

SU x + TU y = 0, 

where 

U = r\S=( c2 - u2 ~ UV \ and T l ~ UV 



1/ \-l 

and put this system in characteristic form formally. This yields the char- 
acteristic equations 

+ A +Xi = 0, K f - A-v^ = 0, 

^ + = 0, - A+Vrj = 0, 

where 

uv ± c ^q 2 - c 2 



A ± = 



c 2 -u 2 



Thus in the subsonic region where q 2 < c 2 the characteristics are com- 
plex. Note that rj are functions of u, v. 

Claim, g, rj may be chosen so that in the real plane the solution is real. 125 

Proof. Let to - u - V— lv, to* = u + V— lv. Suppose £ = g(to, to*) and 
write 

A+(to, to*) = A\(to, to*) + (to, to*) 

A-(cj,co*) = A\(to,to*) - V-Ll2(w, co*), 

with A\ = uv/(c 2 - u 2 ), A2-C ^q 2 - c 2 , Re(q 2 - c 2 ) > 0. If u, v are real, 
then A\, A2 are real and 



A + (to,co*) - A-(to*,co) 



where a bar above denotes the complex conjugate. Hence 

u^(to, to*) = A-(co, co*)vg (to, to*) 
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3. Two dimensional steady flow 



and changing values of the variable to — > to*, to* — > to and because u 
and v are real, we obtain 

u^(co, to*) = A+(co, to*) v^(to, to*) 

or 

u^(to, to*) = A + (to,to*)v^(to,to*). 

This shows n - | is a possible characteristic variable in the subsonic 
region. □ 

A simple example: Consider the Cauchy Riemann equations 

U X + Vy = 0, 
V X - Uy = 0. 

Here A+ - ± V-T and therefore 

+ V^Txf = 0, - "V^Tvf = 0, 

yrj - ^f-lx^ = 0, U-q + V-lVrj = 0. 

So, a> = u - CO - u + are characteristic co-ordinates. 

Note. Various problems could be solved in some limited region, for ex- 
ample, a Cauchy or Goursat problem. Now consider the elliptic re- 
gion. There A.+ are complex but the system remains valid. We consider 
x, y, u, v as complex quantities. Note that q 2 - u 2 + v 2 and c 2 (q) is an 
analytic function of q. Of course, we look for solutions which are real 
for real x, y. 

Remarks. (1) Since the solution is analytic and independent of the 
path the number of actually used real variables can be reduced to 

three. 

(2) A ± are analytic in some slit domain because of the singularities 
at q ± c = 0. This surface q ± c = forms a two dimensional 
manifold in four dimensional space (to, to*). 



3.11. Design by the method of complex characteristics 
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We examine a difference scheme for a Goursat problem and first 
consider the real case. The difference scheme is given by 

p{P) ~ P(Q) + MQ)(q(P) - q{Q)) = 0, 
p(P) - p(R) + A-&MP) ~ 900) = 0, 

where we use p, q as variables. We want to solve for p(P) and q(P), see 127 
Fig. 3.7. 




Fig. 3.7. 



Here the data is prescribed on £ = 0, r/ - (For the Goursat problem 
see Garabedian [ 13 1 pp 1 18-1 19). We can solve for p(P), q(Q) provided 
A + + /I- or equivalently q 2 + c 2 . 

In the complex case, the same argument holds, but note that we 
should take A£ = An in the elliptic region. How do we guarantee that 
we have a real solution x, y for real u, v? Note that the equations for x, y 
are linear and they have real coefficients if u, v are real. Hence Re(x) 
and Re(y) are solutions and these are real. 

In practice, it has proved better to prescribe 

x({,n ) = m+g(n ) 
x(% ,n) = g(rj) + f(g ), 

and to get real solutions in the real plane by choosing /(£) = g(g). 
The following problems arise: 

(i) How to choose / so that a stagnation point appears before a sin- 
gularity? 
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3. Two dimensional steady How 



(ii) How to choose / so that the streamlines of the profile join smooth- 
ly at a trailing edge? 

One usually chooses as a trial f,g the corresponding data for the 
Cauchy-Riemann equations and then adjust. 

Remark . Existence up to a singular point follows by a Cauchy-Kowa- 
lewski type of argument for a Goursat problem. 



3.12 Numerical solution with shocks: Off design 
computations 

We now consider the problem of computing flows for supercritical air- 
foils at off design, e.g., with different velocities at infinity than specified. 
This problem was open for a number of years and the first numerical so- 
lution of a nonlinear mixed equation with shock was given by Murman 
and Cole [32). We treat a similar case. Consider the small disturbance 
equation, (I3.19t . 



which is elliptic for <p x > 0, hyperbolic for (f> x < 0. Suppose we have the 
following data on the boundary of the region (see Figure 3.8): 



<t>x <t>xx + <f>yy ~ 



(i) 

(ii) 

(iii) 



y (x,O)eC~ 

<f>y(X,b) = 

<Px(0,y) and <p x (a,y) are given. 



(3.36) 



The values of (f> y on the shaded segment corresponds to a given shape of 
airfoil y = Y(x) in the small disturbance approximation. 



3.12. Numerical solution with shocks: Off... 



Ill 



<t>x given 



= 



4> v e Co- 
Fig. 3.8. 



4> x given 



Remark . Even for this boundary value problem, there is still no exis- 
tence theorem establishing a weak solution. 

Let Uj represent an approximation for (p at the mesh poit x = iAx, 
y - jAy. The form of difference scheme proposed for the elliptic region 
is: 



1 



U J - U J 

.[(_i±i__l) 2 



( ' : 1-1 ) 2 ] + 



1 



7+1 



2Uj + Uj *] = 0, 



2Ax LV Ax ' v Ax ' A ' (Ay) 2 

(3.37) 

a second order accurate scheme. Note that in this forms the difference 
analogue of j> tp 2 dy - cp y dx - 0. This is in so called conservation form. 
In the hyperbolic region we retard the x-differences and use 



1 



2Ax 



Li- 



Ax 



i '- 1 ) 2 _( 



Ax 



U J 
- 2 ) 2 ] + 



1 



(Ay) 2 



[u%-2ul l + u£i = o 

(3.38) 
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This is accurate to second order for some equation of the form 

\(<&x + <t>yy = «<P%)x (3-39) 

130 where e - for <p x > 0. The scheme thus introduces articicial dissipa- 
tion but only in the supersonic region, cf. Chapter I §§ 5 and 7. 

Exercise . Investigate the one dimensional form of ( 13.391 ) and see how 
the solution depends on the parameter. 

The object is to prescribe appropriate boundary conditions using 
d3.36t and then to solve ( 13.371 ) and d3.38t for U J r For this a particular 
relaxation method is used which we will interpret as a time dependent 
problem on the original equation. Consider 

<Pxt - <f>x<Pxx + <f>yy, 

with the same boundary conditions. Then solve numerically by differ- 
encing the following: 

xt = (f> x (f> xx in nAt ^ t ^ (n + l)At (3.40) 
(j>xt = <j>yy in in + l)Af ^ t ^ (n + 2)At (3.41) 

Here At is small. Define 



in nAt ^ t ^ (n + \)At 

in (n + l)At <t^{n + 2) At 



(3.42) 



This is an alternating direction scheme and it turns out that if the quan- 
tities are smooth, this alternate direction scheme implies 

<f>* -> (p as At — > 0. 

131 However, we do not expect a smooth solution; a shock occurs and across 
this shock <p x increases in accordance with the entropy condition if the 
difference have been retarded in the x - direction and if the shock can 
be described by x = X(y). Furthermore one belives that tp approaches a 
steady state as t — > oo and this steady state is the desired solution. 



3.13. Nozzle flow 
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If we now difference the alternating direction time dependent 
scheme in time and space we get a particular relaxation scheme for find- 
ing Uj. 

Alternative schemes have been suggested by Engquist, Osher [9], 
who show in fact that there exists a solution for their difference scheme. 
The main ideas of this method have been used by Jameson [18] and 
others to solve the full equations and in recent work the method has been 
applied to three dimensional flows where the computing difficulties are 
very great. 

3.13 Nozzle flow 

Another transition flow from subsonic to supersonic occurs in a nozzle 
but with much more stable behaviour. 

The simplest example is a Meyer flow (see Bers [2]) where there is 
an elegant exact solution. From (13.81) one finds the equations for q, 9 as 
functions of (p, tp and one finds near q - c that 

S$ = 9(f„ SS (f, — 9,/,, 

where S is related to cr. For every A, 

A 2 

S =A<p + —it; 2 
A 3 

9 = A 2 (ptp + — l/r 5 
6 

is a solution. The flow in the <p, \p plane has, of course, as streamlines 
the horizontals ip = constant. Thus the wall of the corresponding flow 
in <p, tp plane consists of two straight lines. 

The sonic line is a parabola and the characteristics are (9 - 9 ) 2 = 

4 

-cr 3 and there are four of them passing through (p = 0, tp = given 
A , 

by (p - ±-^r. Hence the mapping into the hodograph plane is not 

one-to-one but there is a fold and the region 9 2 < ^cr 3 is covered three 
times. 
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3. Two dimensional steady flow 



To design a general nozzle with prescribed speed at infinity, we 
therefore solve 

Kl/fgg + tp (T(T = 

in the hodograph plane cr, 9. Specify the singularity at cr^ by analogue 
with incompressible flow (see the Figure 3.9). Leave out the gap made 
by the characteristics through the origin. 




Characteristics 
Fig. 3.9. 




Next continue the flow across the characteristics to get two layers 
ending on the characteristics. Continue the flow on the third sheet into 
the whole quadrant bounded by four characteristics. 

To continue the flow, one uses hyperbolic methods. The flow can 
be terminated by a shock and the outgoing flow from the nozzle will be 
subsonic. 
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